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INTRODUCTION 


Studies in the mathematical modeling of the high-speed turbulent 
combustion has received renewal attention in the recent years. The review of 
fundamentals, approaches and extensive bibliography was presented by Bray, Libbi 
and Williams (1994). 

In order to obtain accurate predictions for turbulent combustible flows, the 
effects of turbulent fluctuations on the chemical source terms should be taken into 
account. The averaging of chemical source terms requires to utilize probability 
density function (PDF) model. There are two main approaches which are dominant 
in high-speed combustion modeling now. In the first approach, PDF form is 
assumed based on intuitia of modelliers (see, for example, Spiegler et.al., 1976; 
Girimaji, 1991 ; Baurle et.al.,1992). The second way is much more elaborate and it is 
based on the solution of evolution equation for PDF. This approach was proposed 
by S.Pope (1981) for incompressible flames. Recently, it was modified for modeling 
of compressible flames in studies of Farschi (1989); Hsu (1991); Hsu, Raji, Norris 
(1994); Eifer, Kollman (1993). But its realization in CFD is extremely expensive in 
computations due to large multidimensionality of PDF evolution equation (Baurle, 
Hsu, Hassan, 1994). 

Promising way for development of computationally non-expensive procedure 
for PDF construction is flamelet approach (FL) formulated for subsonic turbulent 
flames in studies of Bilger (1980), Kuznetsov ( 1977, 1982), Peters (1984) and Williams 
(1975). The simplification of the turbulence/chemistry interaction modeling is 
achieved here based on the assumption that chemical processes are mostly confined 
to the local vicinity of the stoichiometric surfaces. This assumption allows to reduce 
instantaneous conservation equations for reactive scalars to the system of ordinary 
differential equations (flamelet equations). Its solution gives relations for reactive 
species mass fraction and temperature depending on mixture fraction z and its 

scalar dissipation N = D(Vz)~ where D is molecular diffusivity. The FL solutions are 
used to present joint PDF of reactive scalars ,/(C 1 ,...,Cj,T) as a function of mixture 
fraction and scalar dissipation joint PDF ,/(z,N) only. So consideration of the 
reactive scalars statistics is reduced to the modeling of large- and small-scale 
statistics for the passive scalar. Additional benefits of the FL approach are 
connected with the fact that flamelet equations can be integrated before starting 
hydrodynamics calculations for many cases of turbulent flames. So, CFD modeling 
can be performed using tabulated solutions for reactive scalars versus mixture 
fraction z and scalar dissipation N (so-called flamelet library approach). By such a 
manner, very complex detailed kinetics schemes can be incorporated into the CFD 
codes without significant computational time increasing. 

Previously, FL approach demonstrated satisfactory results in predictions of 
reactive species and temperature for subsonic laboratory nonpremixed flames 
(Liew, Bray, Moss, 1984; Buriko, Kuznetsov, 1988; Buriko et.al., 1994). Currently, FL 
model capabilities are studied for simulation of the combustion and pollutants 
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emissions in gas-turbine combustors (Buriko et.al. 1996; Bai, Fuchs 1995; Amin et.al. 
1995; etc.). 

The flamelet modeling of the supersonic flames is much more complex 
compare to the case of low-Mach number flames The main problems which 
complicated FL modeling of compressible flames are connected with: 

(i) Large input of flow kinetic energy term in energy conservation equation and 
pressure variation inside the flow field which provide dependence of the reactive 
scalars on the additional parameters i.e. on flow velocity V and pressure p: 
C a =C tt (z,N,V,p). 

So, dependencies for the reactive scalars became 4-parametric instead of 2- 
parametric ones for low-Mach number flows. This feature complicates procedure for 
flamelet library development and averaging. The dependence of C a on V and p 
worsen the accuracy of hydrodynamics and kinetics splitting also. 

(ii) Large ignition delay. Often supersonic jet flames are not stabilized at nozzle lip 
and ignition takes place in some cross section downstream the nozzle exit. The 
ignition delay distance can be significant for supersonic flames. The basic 
assumptions of the FL approach are violated in the vicinity of self-ignition point 
since reagents can be partially-premixed in self-ignition region and chemical 
processes can occur in some volumes or near-premixed flame fronts. 

Nevertheless, possible benefits, which can provide flamelet approach for the 
simplification of the turbulence/chemistry interaction modeling, requires modeliers 
to look for possible ways for model modifications and its implementation into the 
simulation of combustion processes in supersonic flames. The estimations of 
possible combustion regimes in hydrogen/air systems for application to scramjet 
combustors were reported by Balakrishnan and Williams (1993). Their results 

indicated that, close to flamelet, combustion regimes are likely exist in supersonic 
combustors at typical operational conditions of hypersonic airplanes. 

First studies toward flamelet model generalization for supersonic flames have 
been performed recently by Zheng and Bray (1994), Bezgin et.al. (1995), 

Morgenthaler et.al. (1997). In particular, extension of the Kuznetsov (1982) flamelet 
approach for the expanded supersonic jet flames was done in 1995-1996 under 
implementation of the NASA Cooperative Agreement NCCW-75 by our group 
(Bezgin et.al. 1996). It was found that FL approach allowed to describe satisfactory 
main features of supersonic H 2 /air jet flames. Model demonstrated also high 
capabilities for reduction of the computational expenses in CFD modeling of the 
supersonic flames taking into account detailed oxidation chemistry. However, 

disadvantages and limitations of the existing version of FL approach were found in 
that study also. They were: i) inaccuracy of reactive scalars predictions obtained in 
flamelet modeling for well-known Burrows— Kurkov (1973) test case; ii) significant 
inaccuracy in predictions of the ignition delay distance; iii) improper model 

operation in the vicinity of self-ignition point. 
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So, studies toward further improvement of the FL approach capabilities were 
continued this year under NASA LeRC Cooperative Agreement NCC3-496. They 
were performed in two main directions: 

1. Analysis and improvement of previously found inaccuracy in model predictions 
and consideration of the new test cases for more careful verification of the model 
capabilities. 

2. Searching possible way for modification of the FL approach for its more correct 
operation in the case of flames with large ignition delay distance. 

The current Report summarizes results of the performed studies. It is organized as 
follows: 

Sec. I contains cumulative summary of the flamelet approach which is 
necessary for explanations in the next Sections of the Report. 

Sec. II contains brief description of the test cases which were used both for 
the analysis of the previously obtained discrepancy and in the current CFD tests of 
the flamelet approach. 

Brief analysis of the previous results and current model modifications are 
given in Sec. III. 

The results and analysis of the current CFD tests are presented in Sec. IV. In 
current study, we considered 4 tests cases which covered 3 main types of the flow 
field (round jet, planar wall jet and planar mixing layer) and two classes of flames 
(with and without large ignition delay distance). All the considered test cases were 
calculated based on the averaged system of 2-D Navier- Stokes equations and using 
the same turbulence and detailed chemistry models. The description of the obtained 
results is presented in: 

Sec. IV. 1 - for Burrows and Kurkov (1973) wall-jet test; 

Sec.IV.2 - for Beach round jet test reported in Evans et.al (1978)\ 

Sec.IV.3 - for Cheng et.al. (1994) round jet test; 

Sec.IV.4 - for Chang et.al. (1993) planar mixing layer experiment. 

The supplemented information about CFD tests is submitted to the 
Appendixes (description of the mathematical model, thermodynamics and kinetics 
reference data, etc.). Only brief description of applied computational procedures is 
presented in current Report since we used the same codes as those described in our 
previous Final Report to NASA (Bezgin et.al., 1996). 

Sec.V contains results of our investigations toward modification of the 
flamelet approach for flames with large ignition delay distance. Our reasons and 
modified flamelet equations (MFL) are presented in Sec.V. 1. Computational 
procedure, which was used for their solution, is outlined in Sec.V.2. Results of 
modeling tests of the MFL equations are presented in Sec.V. 3 together with 
suggestions concerning possible way for model improvement. 

The Report is finished by Summary and suggestions for further investigations. 

Research team greatly thanks to our supervisor Dr. Louis Povinelli (NASA 
LeRC) for the formulation of the research directions for current investigation and 
for fruitful discussions. We greatly thanks also to Prof.R.Pitz (Vanderbilt University, 


6 


US), Drs. C. Chang (NASA LeRC), T.Cheng (Chung-Hua Polytechnic Institute, 
Taiwan), O.Jarrett (NASA LaRC) for kindly presented experimental data and for 
fruitful discussions of their test cases. We very thanks also to Prof. S.Pope (Cornell 
University, US) and Prof. V.Sabelnokov (TsAGI, Russia) for fruitful discussions 
during the International Colloquium on Advanced Computation&Analysis of 
Combustion (Moscow, 1997) where part of the current results was reported. 
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I. FLAMELET APPROACH 


The current concept of the flamelet approach (FL) remained the same as that 
presented in Bezgin et.al (1996). The schematic of the approach is given in Fig.l. It 
is based on the assumption that chemical reactions are mostly occur in narrow 
regions (flamelets) located in the vicinity of the surfaces with stoichiometric 
composition. 

The stoichiometric surface is highly curved and randomly fluctuated in 
turbulent flows. So, the mixture fraction z is used to characterize its location. The 
mixture fraction z is introduced as the total mass fraction of all kinds of atoms 
initially been contained in fuel and then converted to other chemical species arising 
in the flame. It obeys conservation equation without source term: 


^ + ^ = V pDVz 

at ax k 


at) 


Value of z = 1 corresponds to pure fuel and z = 0 - to pure air; mixture fraction has 
value z = z s = 1/(1 + St) at the stoichiometric surface (St is the stoichiometric factor). 

The assumption about small thickness of the reaction zones allows to simplify 
modeling of the turbulence/chemistry interaction problem (Fig.l). First of all, it is 
possible to reduce instantaneous conservation equations for reactive scalars to the 
system of ordinary differential equations (flamelet equations). Its solution gives 
relations for reactive species mass fractions and temperature depending on mixture 

fraction z and its scalar dissipation N = D(Vz)" where D is molecular diffusivity i.e. 

C u = C u (z,N), T = T(z,N). These relations are used to present joint PDF for reactive 
scalars i/(C 1 ,...,Cj,T) depending on mixture fraction and its scalar dissipation joint 
PDF .y(z,N) only. Let us present flamelet equations and passive scalar statistics 
model in a little more details. 


Flamelet Equations. Based on assumption about small thickness of the reaction 
zones, the unsteady and convective terms in the instantaneous conservation 
equations for reactive scalars are dropped out and mixture fraction is used as 
independent variable instead of space coordinate. Full account was done by 
Kuznetsov (1982) and it was repeated in our previous Final Report (Bezgin et.al., 
1996). The same kind of reasoning is applied to the energy conservation equation 
also. Resulting system of flamelet equations has the form: 

N s -^-^ + R a =0 (1.2) 

dz" 


dz 2 

where C (t is mass fraction of a-specie (a=l,...,J); J is total number of reactive 
species; R„ is the chemical production term; H is total enthalpy defined as 
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H = h + < VV >; h = ih„ C u is static enthalpy; species specific enthalpies 

2 a 1 

h a = J Cp a dT + Ah a (T 0 ) are used taking into account species heats of formation Ah a 

T 0 

- s ( dzV 

at reference temperature T 0 ; V is velocity vector. The parameter N = D is 

the value of instantaneous scalar dissipation N = D(Vz) 2 at the stoichiometric surface 
which characterizes the reagents fluxes to the reaction zones; n is coordinate 
normal to the stoichiometric surface z = z s (Fig.l). 

The boundary conditions (BC) for the flamelet equations (I.2)-(I.3) are posed 
at z=l (pure fuel) and z = 0 (pure oxidizer): 

z = 0 H = H A ;C a =C* (1-4) 

z=l H = H 1 ;C a = C£ a= 1 J 


where superscripts F and A denote composition and total enthalpies of fuel and 
oxidizer respectively. 

The eq.(1.3) is integrated over z from 0 to 1 and the total enthalpies of fuel 
(at z=l) and air (at z = 0) are used to define the constants in the obtained linear 
relation. As the result, flamelet model equations (1.2), (1.3) are re-written in a form: 

+ R =0 CL— 1.....J (1.5) 

dz 2 

V 2 

h = (H f - H a )z + H a 

2 


The formulated flamelet model boundary problem (1.5) with the boundary 
conditions (1.4) give the solution for C a and static temperature T in the following 
parametric form: 


C u = C„(z, N s , p s , BC); 


a = 1,...,J 


( 1 . 6 ) 


T = T(z, N s , p s , 


X! 

2 


, BC) 


where p s is pressure in the reaction zone and BC denotes boundary conditions (1.4). 
The solution (1.6) is considered in the flamelet approach as an instantaneous 
relations between the reactive species mass fractions and temperature on the one 
hand and mixture fraction, scalar dissipation at the stoichiometric surface, local flow 
velocity and pressure on the other. Additional simplifications were adopted in 
Bezgin et.al (1996) to decrease number of variables in parametric relations (1.6). 
These simplifications were: 


i) the role of the pressure fluctuations on the combustion chemistry was ignored; 
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ii) correlation between flow velocity and mixture fraction distributions was applied 
in eq.(1.5) in a simplified form, which is approximately valid for unconfined jets 
(Abramovich et.al., 1984): 

v-v A p 

y F _ y A ~ z 

where V A ,V F are the mean flow velocities in the air and fuel core flows 
respectively, P is some exponent which was chosen as P«l/Sc t (Sc t is the turbulent 
Schmidt number). 

Recently, practically the same kind of reasoning was applied by Morgcnthaler 
et.al. (1997) in their FL modeling of supersonic jet flames. 

Such simplifications allowed to split flamelet equations from the 
hydrodynamical ones and to apply flamelet library concept. This reduces flamelet 
modeling of reacting flows to the two-step procedure schematically shown in Fig. 2. 
At the first step, flamelet eqs.(1.5) are solved for different values of N s and pressure 
p s . The obtained solutions are collected into the flamelet library in a parametric 
form on z, N s and pressure p s . At the second step, the obtained flamelet library is 
used together with CFD code for the flow hydrodynamics modeling. This two-step 
procedure is described in Appendix A more detailed. 

It is important to note flamelet model behavior in the vicinity of 
ignition/extinction conditions. Let us illustrate it using typical dependence of the 
calculated maximum flame temperature T m vs scalar dissipation N s which is 
schematically shown in Fig. 3. It is seen that increasing of the scalar dissipation N s 
(i.e. growth of reagents mixing rate) lead to increasing of the chemical 
nonequilibriumness in the reaction zone and hence it decreases maximum flame 
temperature T m . It is seen that burning solution exists only if N s < A' ( ". The flame 
quench occurs when N s becomes higher than This means that the flame 

extincts when the mixing rate (which is controlled by scalar dissipation) becomes 
too high compare to reactants consumption in chemical reactions due to the 
limitations of the finite rate chemistry. It is seen also (Fig. 3) that there is the second 
critical value of the scalar dissipation which corresponds to the transition from 
pure mixing to combustion regimes. As a rule, >> Ay; 2 1 . Both N [ ] r 1 and N ' 
values can be obtained from the flamelet calculations for each particular operational 
conditions (pressure, air and fuel temperatures) and given detailed chemistry model 
for fuel oxidation. 

The existence of the critical values of scalar dissipation allows to restrict 
flamelet parametric calculations during flamelet library production by the range 
N s e[0, jV F) ). For the higher values of N s the pure mixing solution can be used to 
calculate the mixture composition and thermodynamics properties. 
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Passive Scalar Statistics Model. The flamelet solutions (1.6) allows to present joint 
PDF for reactive scalars and temperature in a form: 

P(z,N s ,C 1 ,...,C J ,T) = P(z,N s )5(T-T F ')n§(C a -0 

a-1 

where T F1 = T(z, N s ) and C FI =C a (z, N s ) are the solutions of flamelet equations. It is 
seen that only passive scalar joint PDF P(z,N) is needed to obtain mean values of 
reactive scalars. It can be obtained from the numerical solution of the evolved PDF 
equation for passive scalar or using the algebraic models for passive scalar PDF. We 
used the second approach. 

The Favre joint PDF of mixture fraction z and scalar dissipation .^z,N s ) was 
considered in a following form: 

^z,N s ) = ^ P(z, N s ) = (1 - y )5(z)5(N' s ) + yP t (z)5(N s - N? ) 

P 

where y is the intermittency factor; P, is the mixture fraction probability density 
function in a turbulent mixing layer; 5 is the Dirac function. 

The intermittency factor y and probability density function in a turbulent 
mixing layer P t were approximated based on self-similar solutions of PDF equation 
presented in Kuznetsov and Sabelnikov (1990) and used by us previously (Bezgin et 
al. 1996). The intermittency factor y was calculated using approximate relation: 
fl.31/(l + c 2 / z 2 ) if a/z>0.555 

Y = < _ U* > ) 

| 1 if a / z < 0.555 

where z = pz/p is Favre averaged mixture fraction and a : =pz"z"/p is mixture 
fraction variance. The conditional PDF in the turbulent mixing layer P,(z) was taken 
in Airy/Gaussian form reported in Appendix A. 

The conditionally averaged value of scalar dissipation at the stoichiometric 
surface was approximated as: 

_ n| 

Nf _ (1.8, 

n, „ 


where N 


I Z'-ZS 

intermittency 
fraction z = z s 


v| are the mean values of the scalar dissipation N and 
factor y calculated under the condition that mean value of mixture 


The conventional approximation for the mean scalar dissipation N was used: 


where K is turbulence kinetic energy, v t is eddy viscosity, c diss is empirical factor. 

In such a treatment, only turbulence kinetic energy K, eddy viscosity v,, 
mean mixture fraction z = pz/p and its variance cr = pz"z" / p are needed to 
calculate the PDF in any given point of the flowfield. These turbulent mixing 
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characteristics were calculated using conventional semi-empirical transport 
equations of the turbulence modeling (see, Appendix A) which were solved together 
with averaged hydrodynamics equations as it is sketched in Fig. 2. 


Self — ignition criterion. Often supersonic flames are not stabilized at nozzle lip and 
self-ignition takes place in some cross section downstream the nozzle exit as it is 
schematically shown in Fig. 2. The ignition delay distance can be significant for 
supersonic flames and should be modeled. 

The reactants can be partially-premixed in the vicinity of the self-ignition 
point. So, chemical processes can occur, particularly, in some volumes or near- 
premixed flame fronts. Unfortunately, accurate calculations in the vicinity of 
ignition point can not be performed using FL since this model was constructed for 
thin diffusion front (some model modification which can provide more accurate 
modeling in the near-ignition regions is presented in Sec.V of the Report). 
However, upstream (in the mixing region) and downstream (in the diffusion flame 
region) self-ignition region, flamelet approach can be used. That is why, we applied 
somehow simplified approach for prediction of self-ignition point. 

The combustion in the vicinity of the self-ignition point was expected to be 
in diffusion flame mode. Possible existence of the multiple combustion regimes or 
volume exothermic reactions was ignored. The self-ignition point location was 
predicted approximately, as a point on the mean stoichiometric surface, where 
conditionally averaged value of scalar dissipation became lower than its critical 

value N cr . This simplified treatment is schematically shown in Fig. 2. Additional 
explanations regarding applied, in the current study, approach for the prediction of 
the self-ignition point location are presented in Sec. III. 
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II. SUMMARY OF THE TEST CASES 


Complexity of the measurements in compressible flames significantly 
decreases number of the test cases which can be used for CFD validation of the 
turbulent combustion model. The analysis of the known unclassified literature and 
consultations with our supervisor Dr. Povinelli (LeRC) allowed to select 4 test cases 
where sufficiently detailed information about test conditions, flowfield and reactive 
scalars characteristics was reported. They cover 3 main types of the flow field 
(round jet, planar wall jet and planar mixing layer) and two classes of flames (with 
and without large ignition delay distance). The brief summary of the selected test 
cases is presented in Table 1. More detailed information and sketches are presented 
in Sec. IV which reports final results of current CFD. 


Table 1. TEST CASES SUMMARY 


SOURCE 

CASE 

FLAME 

Mh2 

M a 

Th2 

[K] 

T air 

[K] 

Pm/Pair 

[MPa/MPa] 

Ignition 

delay 

TOOLS 

Burrows & 

Kurkov 

(1973) 

NASA 

TM 

X-2828 

planar 
wall -jet 

H2-air‘> 

mixing 

1 . 

2.44 

254 

1150 

0. 1/0.1 


P' P 

GS, TC 

H 2 -air** 

flame 

1 . 

2.44 

254 

1270 

0. 1/0.1 

large 

Beach et.al. 
(1978) 

NASA 
TP 1 169 

round- 

jet 

H 2 -air # * 

flame 

2. 

1.9 

251 

1495 

0. 1/0.1 

small or 
none 

p\ 

GS, 

Cheng et.al. 
(1994) 

Comb&Fl. 

v.99 

round- 

jet 

H 2 -air‘> 

flame 

■ 

2. 

545 

1250 

0. 1/0.1 

mid 

Spont. 

Raman 

L1F 

Chang et.al. 
(1993) 

A1AA- 

93-2381 

planar 

mixing 

layer 

air-air 

mixing 

0.39 

0.72 

303 

824 

0. 1/0.1 


LDV 

TC 

h 2 /n 2 - 

air 

flame”* 

0.3 

0.71 

348 

817 

0.106/0.106 

pilot 

flame 


’> vitiated air (obtained by burning of H 2 in enriched by 0 2 air) 

") pilot H '2 flame in air flow for combustion process stabilization 

Tools are diagnostics techniques: P r - pitot pressure probes, P w - wall static pressure measurements; 
GS-gas sampling, TC - thermocouple. 
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III. PREVIOUS TESTS and CURRENT MODEL MODIFICATIONS 

The Beach experiment reported in Evans et.al. (1978) and Burrows, Kurkov 
(1973) test were considered in previous studies (Bezgin et.al. ,1996). The mostly 
important, for further explanations, results of the previous tests are summarized in 
Fig. 4. Here FL predictions for water concentration are plotted together with data for 
both considered previously test cases. It is seen that accuracy of FL predictions was 
quite different. The FL predictions were quite satisfactory for Beach test case. 
However, unsatisfactory results were obtained for the Burrows-Kurkov experiment. 
The significant discrepancy (about 40%) between experimental data and results of 
computations was found not only in H 2 0 concentration peak value but in static 
temperature distributions and ignition delay distance (predicted 0.08 m downstream 
fuel injector instead of *0.18m reported by Burrows and Kurkov). Trying to improve 
accuracy of predictions, we performed a large number of draft parametric tests, 
where we varied turbulence models (Secundov's one parametric "v t -90 model 
presented in Appendix A, standard k-e model of Jones&Launder (1973) and Chein 
(1982) k-e model were used). We varied also shape of the inflow profiles and level 
of the turbulence intensity at the entrance of the test section. Unfortunately, all 
these modifications did not provide any reasonable improvement of the predictions 
(or even made them worse). Nevertheless, analysis started in (Bezgin et.al., 1996) 
allowed to clarify two possible roots of the discrepancy. They are: (i) too high level 
of the mixture fraction fluctuations intensity INT = a/ z predicted by the turbulence 
model equation for mixture fraction variance and/or (ii) too high chemical 
nonequilibriumness predicted by the used detailed kinetics model. The relative role 
of these two effects on the previous (and unsatisfactory) model predictions for 
Burrows&Kurkov experiment is illustrated by Fig. 5 taken from Bezgin et.al. (1996). 
Here, the measured by Burrows and Kurkov, H 2 0 mole fractions (denoted by 
crosses) are plotted versus mixture fraction z. The dashed line denotes the 
equilibrium chemistry solution for water mole fraction x‘h:o plotted vs mixture 
fraction z also. It was calculated based on the assumption that both chemical 
nonequilibriumness and scalar field fluctuations are absent in the flow so it is the 
upper limit for possible water concentration distributions. It is seen that 
experimentally measured H 2 0 mole fraction distribution is very close to the 
equilibrium limit XnVo and so role of both chemica l nonequilibriumness and 
fluctuations intensity were not too high in this conditions. The flamelet model 
overpredicted role of these effects as it is shown in Fig. 5. The instantaneous 
flamelet model solution obtained in calculations for water mole fraction 
XH20 = XH2o(AN s ,p s ) is plotted vs mixture fraction z by fine solid line. The averaged 
distribution of water mole fraction: 

X„», =Jx„ ! „^,N s ) dzdN s 

is plotted by fat solid line. It is seen approximately 50%:50% input of both chemistry 
nonequilibriumness and averaging procedure into the resulting averaged water 
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mole fraction distribution Xh 2 o predicted by the FL approach. It is seen that 
averaged concentration peak value is 40% lower than the peak value for X 1 I 20 
distribution. Such a predictions were very far from the experimental data of Burrows 
and Kurkov. 

Previously we connected this model uncertainties with improper operation of 
only mixture fraction variance equation in our turbulence model. However, 
additional tests of the detailed chemistry approximations performed in the current 
study required us to change this point of view and to introduce modifications not 
only into the mixture fraction variance equation but into the detailed chemistry 
model also. Trying to improve flamelet model predictions for self-ignition point 
location, we changed self-ignition criterion also. Our reasoning and account of 
introduced modifications are presented below. 


Mixture Fraction Variance Equation. To decrease intensity of the mixture fraction 
fluctuations we somewhat changed values of empirical coefficients in balance 
equation for mixture fraction variance a 2 =pz"z"/p. Its general form remained the 
same as that used by us previously and widely distributed between modellers (see, 
for example, Jones&Whitelaw 1982): 


dp a 2 
at 


+ VpVa 2 


a 

dx,. 


P c a v t + 


Sc 



diffusion 


+^~M Vz ) 2 - 2 p n 

production dissipation 


(in i) 


where c a is empirical constant which is usually chosen as c CT «l/Sc t ; N is mean value 
of the scalar dissipation. 

The changes were made in the value of empirical factor in approximation for the 
mean scalar dissipation. We used in our previous predictions the following relation: 


N = c 


diss 


K-ct 2 

v t 


(HI 2) 


where K is turbulence kinetic energy, v t is eddy viscosity, c djss is empirical factor. 
Previously, we used value of c dtss = 0.07. 

Somewhat different form for approximation of the mean scalar dissipation is 
commonly used by modellers: 


K 

where e is the turbulence energy dissipation rate and value of empirical factor c is 
usually chosen as c*0.95-l (Jones&Whitelaw, 1982; Bilger, 1989; etc.). So, one can 
obtain that c djss «0. lc based on conventional, in turbulence modeling, relation 
K 2 

e w 0.1 . That is why, we somehow increased role of the dissipative term (and 

v t 
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decreased intensity of fluctuations) by increasing coefficient c djsK in (III. 2) from 
previously used value 0.07 to 0.1. 


Detailed Chemistry Approximation. Previous tests of the flamelet model were done 
using slightly simplified hydrogen oxidation block taken from the detailed kinetics 
scheme proposed by Miller and Bowman (1989). The resulting detailed kinetics 
model included 21 reactions between 11 species (H 2 , 0 2 , H 2 0, H, O, OH, H 2 0 2 , 
H0 2 , N 2 , N, NO). The used, previously, simplification was that the difference in 
third-body efficiencies for three-molecular reactions (like H + H + M = H 2 + M, 
OH + H + M = H 2 0 + M, etc.) were ignored. It was expected that efficiencies of all 
molecules, which act as third body M in three-molecular reactions, are the same as 
that of molecular nitrogen N 2 . Previously, we expected that such simplification did 
not significantly influence results of chemistry calculations. Unfortunately, 
additional tests showed that such simplification of the detailed chemistry model is 
too rough. 

We performed separate flamelet model calculations for typical values of the 
scalar dissipation in the reaction zone which corresponded to mid and far field of 
supersonic jet flames (N s = 5-50sec''). Three detailed kinetics mechanisms were used 
in sensitivity tests. They were: 

(i) previously used simplified H 2 oxidation block from Miller&Bowman 

kinetics scheme; 

(ii) the same H 2 oxidation block from Miller&Bowman scheme with different 

third body efficiencies which were taken into account in accordance 
with the data reported by Miller and Bowman. 

(iii) Jachimowski (1988) scheme, which was developed in NASA for 

modeling of H 2 /air oxidation chemistry in supersonic flames. This 
model utilizes different third body efficiencies also. 

The kinetic schemes and third-body efficiencies are presented in Appendix C. 

Typical results of flamelet calculations for conditions of Burrows- Kurkov 
experiment are presented in Fig. 6. It was found that introduction of the different 
third body efficiencies should increase equilibriumness of the flame in the regions 
of "well-developed" combustion. It is seen that Miller&Bowman scheme with 
selective third-bodies efficiencies predicted somehow higher flame 
nonequilibriumness compare to the Jachimowski scheme but, in general, Miller- 
Bowman and Jachimowski detailed kinetics schemes gave close to each other 
results. 

It was found also that such detalization of the kinetics scheme significantly 
influenced predictions of critical values of the scalar dissipation N cr . The summary 
of the and values obtained in sensitivity tests is presented in Table 2. 
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Table 2. DETAILED KINETICS TESTS 

(Burrows&Kurkov test case conditions) 


Kinetics model 

N (U 

sec' 1 

n (2) 

1 cr ' 

sec" 1 

Miller&Bowman scheme 
(selective efficiencies are ignored) 

641.0 

89.6 

Miller&Bowman scheme 
(selective efficiencies are used) 

942.9 

14.2 ! 

Jachimowski scheme 
(selective efficiencies are used) 

1091.4 

15.2 


The same influence of the third bodies selective efficiencies was obtained for 
other test cases also. 

The obtained result required us to introduce selective third body efficiencies 
into the detailed chemistry model. Additionally, we choose Jachimowski H 2 /air 
oxidation scheme for final CFD tests since it was specially developed for supersonic 
flames modeling (it was tested by Casimir Jachimowski (1988) using data for high 
enthalpy flows at the conditions close to the regimes which were investigated in the 
current study). 


Self— Ignition Criterion. The previously found (Bezgin et.al., 1996) significant 
underprediction of the ignition delay distance required us to modify self-ignition 
criterion. 

As it was reported in Sec. I, the self-ignition point location is predicted in 
flamelet approach approximately, as a point on the stoichiometric surface, where 
scalar dissipation N * became lower than its critical value N cr ( A,‘ s <N cr ). However, 
the flamelet solutions provide two values of critical scalar dissipation (AT' J and 
jV'; } ) as it was discussed in Sec.I and schematically shown in Fig. 3. The value of 
Ah)- 1 corresponds to the flame quenching and corresponds to the 

mixing/burning transition. As a rule as it is illustrated by data 

presented in Table. 2. So, using N ( c l r ’ or N™ as an indicator in self-ignition criterion 
will lead to different results regarding ignition point location. 

Previously, we used value of as an indicator for ignition point location. 
However, additional analysis led us to the assumption that it is much more correctly 
to use value of N ( J } for prediction of self-ignition point location in supersonic 
delayed flames. Our reasoning regarding this choice is presented below. 
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The using of as an indicator of critical conditions is widely distributed in 
flamelet modeling (Liew et.al, 1984; Morgenthaler et.al.,1997; etc.). The existence of 
the second critical value of scalar dissipation is usually ignored. Such a choice is 
justified for low-speed flames at near-room or slightly elevated conditions. As a 
rule, at such conditions, combustion wave initialization reguires external ignitor 
(spark device or pilot flame). Flameholders are usually applied for combustion wave 
stabilization. That is why, flame is initiated and stabilized by external tools and so 
combustion model should control possible flame quenching due to too high mixing 
rate only. The using of is justified in such situation. Moreover, if flame ignites 
in some point of the low-speed mixing layer far downstream fuel injector, it can 
move upstream up to the fuel injector or to the point where conditions for its 
extinction will occur (so-called "flashback" effect). As the result, the steady-state 
combustion zone would exist in the whole region downstream the point where 
scalar dissipation became lower than N ( J r } . 

One can expect that this reasoning is violated for self-ignition in supersonic 
flows since flame movement upstream self-ignition point is quite questionable due 
to very high value of the flow velocity. So, we expected that it is more justified to 
use for approximate prediction of self-ignition point location in high-speed 

flows since it is the parameter N ( J ] responsible for mixing/burning transition of the 
flamelet solutions. 

Final CFD tests (Sec. IV) showed, that using of in self-ignition criterion 

provided much better predictions. However, first attempts to apply N'j’ as an 
indicator in self-ignition criteria were unsatisfactory. The reported in Table 2 values 
of lead to significant overprediction of the ignition delay distance for the 

Burrows-Kurkov test case. The same results were obtained for Cheng et.al. (1994) 
test case conditions also. This discrepancy required us to introduce one more 
modification into our flamelet modeling of supersonic flames which is described 
below. 


Vitiated Air Composition. The important feature of three, among considered test 
cases, is the air flow pre-heating by the burning of the H 2 in air enriched by 0 2 . As 
a rule, this leads to the presence of non-zero radicals concentrations (H,O p OH) in 
the high-temperature vitiated air flowing into the test section. Previously , we did 
not take into account this feature. We considered vitiated air flow consisted of 
stable species only. Our reasoning was based on the results of FL calculations 
which showed that both burning solutions of FL equations for reactive species and 
scalar dissipation at quench limit were practically insensitive to such 

modification of air composition. However, this feature influenced parameter N [~ 1 
and hence prediction for self-ignition point location as it is illustrated by Table 3. 
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Table 3. VITIATED AIR COMPOSITION INFLUENCE 

(Jachimowski scheme, Burrows&Kurkov test case conditions) 


Vitiated air composition 

N (l) 

sec' 1 

n (2) 

1 cr ' 

sec" 1 

Stable species only 

(0.258 O 2 + 0.486N 2 + 0.256H 2 O by mass) 

1091.4 

15.2 

Equilibruim composition’) 
(0.258 02 + 0.486N2 + 0.256H 2 0 + 
+ 10 5 OH + 2.8- 10" 4 NO by mass) 

1089.7 

33.8 


'•Species with mass fraction <10‘ 5 arc not shown 


Unfortunately, data about radicals concentrations were presented in only one 
(Cheng et.ai, 1994) among the considered tests. So, to model this effect in the other 
test cases, we considered vitiated air flow as mixture with equilibrium composition. 
It did not change (with accuracy about 0.05%) concentrations of the main species 
compared to the values reported by authors for the considered experiments but it 
led to the presence of small amount of radicals and nitric oxides (with mass fraction 
about 10~ 4 -10' 5 ) in the vitiated air flow. We expect that it is more realistic approach 
than modeling of the vitiated air flow as a mixture of stable species only. 


In general, input of each among the introduced modifications on the results of 
predictions was not very high but their collective influence allowed to obtain 
reasonable improvement of the flamclct predictions in the current tests (Sec. IV) 
compare to those reported previously in Bezgin et.ai. ( 1996). 
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IV. CURRENT RESULTS OF CFD TESTS 


CFD tests of the flamelet approach were done for four test cases summarized 
in Sec. II. The same as previously used, flamelet library concept sketched in Fig. 2 
was applied in current CFD. The modeling of each test case included two 
sequential steps. At the first step, flamelet equations (1.5) were solved parametrically 
and obtained solutions were collected into the flamelet library. At the second step, 
flamelet library was used together with the CFD solver for the flowfield calculations. 
All current CFD tests were done using full system of averaged 2D Navier-Stokes 
equations. The Secundov's one-equation model "V(-90" was used for turbulence 
modeling. Full description of the mathematical model and coupling procedure is 
presented in Appendix A. 

In flamelet calculations, hydrogen combustion chemistry was approximated 
by the detailed H 2 /air kinetics scheme proposed by C.Jachimowski (1988). The 
Jachimowski detailed kinetics model (please, see Appendix C) includes 33 reactions 
between 13 species (H 2 , 0 2 , F1 2 0, H, O, OH, H0 2 , H 2 0 2 , N 2 , N, NO, N0 2 , HNO). 
The species specific enthalpies were used in a form of 7-th order polynomial 
approximations over temperature taken from Alemasov et.al. (1971) and presented 
in Appendix C also. The flamelet solver FLSLV and all details of the flamelet 
libraries generation procedure were fully the same as those reported in Bezgin et.al. 
( 1995, 1996). As previously, flamelet libraries were verified by the interpolation 
accuracy and Richardson extrapolation to zero-length grid step (both were better 
than 1%). The summary CPU requirements for flamelet libraries production and its 
accuracy tests did not exceed 5-10% of total CPU required for CFD modeling. 

The CFD modeling was done using the same modified version of FNAS2D 
code which was used by our group in previous tests of the flamelet approach. It was 
described in our previous Report (Bezgin et.al. 1996). 

The same, as applied in Bezgin et.al., (1996) smoothing of the heat, released 
in the vicinity of the self-ignition point was used in current calculations to avoid 
nonphysical flow disturbances of the flowfield in these regions and to provide better 
convergence of computations. 

The main part of CFD tests was done using HP 9000/735 workstation except 
modeling of the Chang et.al (1993) test case which was modeled using conventional 
PC Pentium 200MHz. 
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IV. 1 BURRO WS-KURKOV WALL-JET EXPERIMENT 


IV. 1.1. TEST CASE 

Scheme of the setup and test conditions of Burrows and Kurkov (1973) 
experiment are given in Fig. 7. The test section was rectangular duct having the 
constant width (0.051m). The air supply duct had 0.089m height. Hydrogen was 
injected parallel to the vitiated air flow. It was injected with a sonic speed through 
the two-dimensional slot located at the backward step in the initial cross section. 
Slot height was h = 0.004m. Lip thickness at the top of the step was 0.76- 10~ 3 m. Test 
section total height expanded linearly from 0.0938m in the initial cross section to 
0.105m at the exit cross section. Composition measurements were done at the exit 
plane of the test section located at x = 0.356m downstream the injector location. 

As it was discussed in Sec. III., this test case was studied previously and we 
met serious difficulties in flamelet modeling. So, the goal of the current calculations 
was both to study how the introduced modifications changed results of predictions 
and to get final conclusions about model accuracy. 

IV. 1.2. COMPUTATIONS 

All details of the computational procedure were not differ from that used in 
our previous modeling of Burrows&Kurkov test case presented in Bezgin et.al. 
(1996). 

The flow field was expected 2-D. The role of the boundary layer at the upper 
duct wall was neglected The simplest molecular transport model was applied i.e. 
the mixture molecular viscosity and diffusivity were estimated based on H 2 
molecular diffusivity and fixed laminar Prandtl and Schmidt numbers Pr = Sc = 0.72. 
The values of the turbulent Schmidt number and constant c„ in turbulent diffusivity 
term of mixture fraction variance eguation (III. 1 ) were Sc t = 1 and c o = 0.67. 

The computational domain together with summary of the boundary 
conditions is given in Fig. 8. The left boundary was located in the hydrogen 
injection cross-section where all parameters distributions were supposed to be 
known in all opened parts of cross-section excluding the slot lip. No-slip velocity 
conditions were posed on the lip wall and on the lower wall of the duct. It was 
expected also that the lip surface and lower wall were adiabatic with zero 
temperature gradient. Hie turbulent kinetic energy and turbulent viscosity were 
equal to zero at the lower wall and lip. The mean mixture fraction and mixture 
fraction variance normal derivatives were equal zero at the lower wall and lip. The 
upper wall was considered as inviscid with zero transversal velocity component. All 
normal derivatives, which were needed to estimate viscous stresses and 
corresponding diffusion fluxes on the wall, were equal to zero. The so-called drift 
boundary conditions with normal derivatives of all parameters determination from 
computational domain were posed in the exit plane of computational domain. 
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The parameters distributions at the inlet boundary were obtained by the same 
manner as it was done previously in (Bezgin et.al., 1996). The longitudinal velocity 
and eddy viscosity distributions in the incoming air flow were derived in 
accordance with the procedure described in Appendix B and using estimated from 
the experimental data boundary layer thickness (8«0.016m). The same procedure 
was used for calculations of u, v t and K distributions in the exit plane of the 
hydrogen injection slot. The 2% velocity fluctuations were expected in free air 
stream. The non-dimensional initial distributions of the longitudinal component of 


the flow velocity U = 


— ; eddy viscosity 

K 



v — 

— — ; turbulent kinetic energy K = 

u a h 


K 


and total enthalpy H = — are given in Fig. 9 where u a is the flow velocity in the air 

K 

core flow and h is hydrogen injector slot height. The uniform step initial profile was 
expected for mixture fraction z. The transverse component of the velocity v and 
mixture fraction variance a 2 were expected to be zero at the entry boundary. 

The reliability of the adopted initial distributions was studied previously (see 
Bezgin et.al. 1996) for non-reacting counterpart of Burrows-Kurkov experiment. The 
satisfactory agreement was found. 

Calculations were done using the same nonuniform grid with 100 cells in 
longitudinal direction and 90 cells in transversal direction (Fig. 8). Grid was 
clustered to the lower wall and to backward-facing step in accordance with 
geometrical progressions. Detailed description of the grid generation procedure was 
given in our previous Final Report also. Both previous and current tests of the 
discretization influence on the results of calculations showed that it was confined in 
small (about 5%) variation of the peak values of turbulence characteristics. 'Idle flow 
hydrodynamics parameters were practically insensitive to the grid variation. 

The convergence was estimated by the L 2 norm for the residual of continuity 
eguation. The L 2 norm behavior vs iteration number did not changed compare to 
the previous results reported in Bezgin et.al. ( 1996). It decreased by 4 orders during 
600 iterations which required about 2 hours CPU of workstation HP 9000/735. 


IV. 1.3 RESULTS and DISCUSSION 


Ignition delay distance. The obtained Mach number field is given in Fig. 10. It is 
seen that the flamelet model predicted self-ignition point location »0.23m 
downstream hydrogen injector. This prediction correlates with accuracy about 20% 
with data reported by Burrows and Kurkov (x«0.18m). Additional tests showed high 
sensitivity of the self-ignition point location to the value of radical concentrations in 
the vitiated air flow. The OH concentration increasing from eguilibrium value 10' 5 
to 10‘ 4 allowed to displace self-ignition point practically into the same position 
where it was observed by Burrows and Kurkov. However, information about radicals 
concentration in vitiated air flow was absent in specification presented by Burrows 
and Kurkov. Their possible level would depend on both pre-heater and air nozzle 
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constructive features. So, we decide to use in our further predictions equilibrium 
values of radicals concentrations in vitiated air as the lower but justified limit. 

Previously, flamelet model prediction was x = 0.08m i.e. the ignition delay 
distance was twice underpredicted. The modifications of the detailed chemistry 
model and ignition criterion improved correlation for the location of self-ignition 
point. It was predicted with the accuracy about 20% (at least). 

Reactive scalars and temperature. The comparison of the obtained reactive species 
mole fraction distributions (solid lines) with the experimental data is given in Fig. 11 
for test section x = 0.356m where the composition measurements were done by 
Burrows and Kurkov. The results of our previous predictions reported in Bezgin 
et.al. (1996) are plotted by dashed lines. It is seen that introduced modifications 
improved model predictions of the reactive scalars. The same improvement was 
obtained for the static temperature distribution also (Fig. 12). The peak values of 
both H 2 0 and static temperature were correlated with the accuracy about 10% 
instead of 35-40% discrepancy obtained previously (Bezgin et.al., 1996). 

The additional analysis showed that predictions improvement was obtained 
by the collective influence of all modifications introduced into the current version 
of the approach. The relative input of the chemical nonequilibriumness and mixture 
fraction fluctuations provided by the current version is illustrated in Fig. 13. The 
relative influence of these factors obtained in the previous unsatisfactory tests is 
presented in Fig. 13 for comparison also. It is seen that both decreasing of the 
fluctuations intensity (approximately by 50% in near stoichiometric zones as it is 
reported in Fig. 14) and increasing of the chemistry equilibriumness provided 
approximately 50% by 50% improvement of the results for reactive scalars. 

It is necessary to note that decreasing of the mixture fractions fluctuations 
intensity was obtained not only due to changes in empirical coefficients of the 
turbulence model equation for mixture fraction variance but due to improvement in 
positioning of self-ignition point also. Self-ignition in conditions of Burrows-Kurkov 
experiment was accompanied by essential transversal expansion of the mixing layer 
in the vicinity of ignition point due to heat release as it is shown in Fig. 15a. Mixing 
layer became much thicker downstream the ignition point and, as a result, mixture 
fraction transversal gradients became lower. Some additional turbulence was 
generated in the vicinity of the ignition point (Fig. 15c). These effects lead to the 
decreasing of the mixture fraction fluctuations intensity (Fig. 15b) through the 
decreasing of the role of production term (which depend on the mean mixture 
fraction gradient) and increasing of the dissipation term in the equation for mixture 
fraction variance. 

We consider the obtained results as quite encouraging however even the 
reported version of the model somehow underpredicted (about 10%) the mainly 
important indicators of heat release, such as H 2 0 and temperature peak values. So, 
we tried to obtain better correlation. Unfortunately all these attempts were 
unsuccessful. We expect now that further improvement of predictions for the used 
version of the flamelet model was prohibited due to very large ignition delay 
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distance leading to possible input of the partially-premixed combustion mechanisms 
(volume reactions or multiple flame configurations) additionally to the main 
diffusion flame mode. These our assumption is based on the comparison of 
calculated and measured wall pressure distributions. 

Wall pressure distribution. Additionally to the gas-sampling and temperature 
measurements, Burrows and Kurkov measured pressure distributions along the 
lower wall of the duct. Unfortunately, as a rule, modeliers do not present results 
provided by their models for this distribution. We know the only publication of 
Kolesnikov (1981) who presented such predictions provided by his assumed PDF 
modeling with Spiegler et.al. (1976) averaging procedure. 

We compared wall pressure distribution provided by current flamelet 
approach with the data of Burrows and Kurkov. It is presented in Fig. 16. 
Experimental data variation in x<0.1m region was attributed by Burrows and 
Kurkov to some uncertainties in positioning of the test section elements and 
existence of the shock wave from the air nozzle so it could not be modeled in our 
calculations. However, it is seen that monotonic pressure increasing downstream 
self-ignition point was obtained along the wall in experiments of Burrow and 
Kurkov. The flamelet model gave rapid pressure rise in the vicinity of the self- 
ignition point (x»0.2m) with the pressure decreasing in the far field. It indicates that 
model overpredicted heat release in the vicinity of the self-ignition point and 
underpredicted it in the far field. Such difference can be explained by improper 
model operation in the vicinity of the self-ignition point. Flamelet predicts very 
rapid transition from mixing to diffusion flame regime and so, it predicts that 
partially premixed reagents should rapidly react in the nearest vicinity of the self- 
ignition point. At the same time it is known that their consumption can realize in 
premixed flames which formed in fuel rich and fuel lean zones additionally to the 
main diffusion front (so-called multiple flame configurations studied ,in particular, 
by Kioni, Rogg, Bray and Linan, 1993) and/or due to volume exothermic reactions. 
The diffusion flame mode input into the total heat release rate would become 
dominant in such a case only sufficiently far downstream self-ignition point. 

We tried to estimate region where such element of partially-premixed 
combustion could be important for conditions of Burrows-Kurkov experiment. For 
this purpose, we varied in our calculations length of domain where we redistributed 
heat initially released by flamelet model in the vicinity of self-ignition point. The 
calculated wall pressure distribution became close to the experimental one only 
when we re-distributed this heat along the whole combustion region up to the exit 
cross-section (see dashed line in Fig. 16). This indicates that the possible input of 
the partially premixed combustion mechanisms can be important up to the exit of 
the test section. 

Based on results of current CFD, it was concluded that introduced 
modifications allowed to improve flamelet approach predictions significantly. 
Reactive species and temperature distributions were correlated with the accuracy 
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about 10-15%. We obtained also much better prediction of self-ignition point 
location. However, additional indirect estimations showed that role of the partially 
premixed combustion mechanisms could be important for this test case and further 
modification of the flamelet approach is needed for more accurate operation in such 
cases. We expect that this feature worsen correlation between reactive scalars and 
mixture fraction and required us to use slightly different values of empirical 
constants (Sc t and l/c„) in diffusivity terms of turbulence model equations for mean 
and variance mixture fraction equations. 

Additional estimations for conditions of Burrows -Kurkov test case were done 
using modified FL equations and they are presented in Sec.V of the Report. These 
estimations correlated with the formulated conclusions. 
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TV.2 BEACH ROUND-JET EXPERIMENT 


I V. 2. 1. TEST CASE 

The sketch of the Beach test case reported in Evans et.al.(1978) is given in 
Fig. 17 together with the inflow conditions. The hydrogen was injected through 
supersonic axisymmetric nozzle with the Mach number M H2 = 2. The hot air was 
obtained by burning of hydrogen in air enriched by oxygen. High-enthalpy vitiated 
air flow was expanded through supersonic nozzle with the Mach number M a =1.9. 
The hydrogen injector tube had external diameter dj = 0.009525m with a lip 
thickness 0.0015m. The air nozzle free stream diameter D was 0.0653m. 

This test case was investigated in our previous study where CFD modeling 
was done using PNS approximation of the hydrodynamics equations. The Beach 
experiment did not meet any serious difficulties in our previous tests except some 
underprediction of the H 2 0 mass fraction for the nearest to the injector test section 
(x/dj = 8.26). However, modifications introduced into the approach required us to 
re-calculate it in the current tests also. 

IV.2.2. COMPUTATIONS 

The assumption about H 2 jet in co-flowing infinite air stream was adopted in 
CFD modeling. The role of the vitiated air flow mixing with ambient air was 
neglected. The flowfield was expected axisymmetrical. The mixture molecular 
viscosity and diffusivity were estimated based on H 2 molecular diffusivity and fixed 
laminar Prandtl and Schmidt numbers Pr = Sc = 0.72. The calculations were done 
for the value of the turbulent Schmidt number Sc t = 0.7 and the value of constant 
c CT = 1/Sc t in eq.(III.l) for mixture fraction variance. 

The schematic of computational domain is given in Fig. 18. The left boundary 
was located in the hydrogen injection cross-section where all parameters 
distributions were supposed to be known in all opened parts of cross-section 
excluding the injector lip. No-slip velocity conditions were posed on the lip wall. It 
was expected also that the lip wall was adiabatic and non-catalytic. The turbulent 
kinetic energy and turbulent viscosity were equal to zero at the lip. The mean 
mixture fraction and mixture fraction variance normal derivatives were equal zero 
at the lip also. The upper boundary of computational domain was at y/dj =2 
position. The no-reflection conditions were posed at the upper boundary. The 
symmetry conditions were posed at the axis of symmetry. The drift boundary 
conditions with normal derivatives of all parameters determination from 
computational domain were posed in the exit plane of computational domain. 

The calculations were done in two sequential domains sketched in Fig. 18 to 
increase computations accuracy in the vicinity of the injector. The first 
computational domain was located in 0<x/djM.5 region and the second one was in 
1.5<x/dj<27.9 region. The drift boundary conditions were applied at the exit plane 
of the first computational domain. All parameters on the inflow boundary of the 
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second computational domain were taken from the calculations in the first domain. 
The computational grid in the first domain was clustered to the lip in both 
transversal and streamwise directions. The grid nodes were clustered to the mixing 
layer in both computational domains. 

The parameters distributions at the inlet boundary of the Domain I were 
obtained by the same manner as it was done in the case of Burrows-Kurkov test 
case. The longitudinal velocity and eddy viscosity distributions in the incoming air 
flow were derived from the estimations of the boundary layer thickness (5 = 0.002m) 
in accordance with the procedure described in Appendix B. The same procedure 
was used for calculations of u, v, and K distributions in the exit plane of the 
hydrogen injection slot. The 2% velocity fluctuations were expected in free stream 
outside the boundary layers. The uniform step initial profile was expected for 
mixture fraction z. The transverse component of the velocity v and mixture fraction 
variance a 2 were expected to be zero at the entry boundary. 

The computations with typical 100x90 grids in each computational domain 
required about 1 hour of HP9000/735 work station. The L 2 norm behavior for the 
residual of continuity equation vs iteration number is presented in Fig. 19. 

The same as in the case of Burrows-Kurkov tests case, series of calculations 
with separate increasing of nodes number in longitudinal or in transversal 
directions was done to control accuracy of computations. In the first series, the 
main computational domain (Domain II in Fig. 18) was divided into 4 subregions 
and calculations were done with patched grids consisted of 100x90 nodes in each 
sub-region. So, the total number of the grid nodes in longitudinal direction was 
increased from 100 to 400 to estimate influence of the discretization in longitudinal 
direction. Additional series of calculations was done using 200x170 (for Domain I 
and 100x170 (for Domain II) adaptive grids. The grid adaptation was realized in 
accordance with spring analogy reported by Baruzzi (1993). In each cross-section all 
grid nodes were supposed to be connected by springs with the stiffness 
proportional to the gradient of Mach number. The obtained adaptive grids are 
presented in Fig. 20 (for Domain I) and in Fig. 21 (for Domain II). 

These tests showed small influence of discretization in both longitudinal and 
transversal directions on the results of calculations. It mostly confined in 2-3% 
variation of the turbulence energy and mixture fraction variance peak values. All 
final results reported below were obtained using presented adaptive grids. 

IV.2.3. RESULTS and DISCUSSION 

Flowfield near injector. Our previous CFD modeling of the Beach test case was 
done using PNS approximation of the hydrodynamics equations. So the flowfield 
features in the vicinity of the injector lip were ignored and PNS marching 
calculations were started from x/dj — 0.33 cross section with some presumed inflow 
profiles of parameters. In the current study we performed CFD modeling using full 
system of averaged Navier-Stokes equations. So, we were able to calculate flowfield 
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in the vicinity of the injector and to control reliability of the inflow conditions 

ad ° Pt The" TtataedTo^rtSd patterns (Mach number field and static pressure 
contours) in the domain near injector are presented in Fig.22. It is se 
compression and expansion waves which arise due to turning and mterachon 
supersonic flows near injector lip. The parameters distributions in x/d,OM^ 
action obtained in current NS computations, were close to presumed 
distributions which were used as inflow conditions for our previous PNS mode mg. 
So detalization of the computational model for near-injector region did not 
significantly influence results of the model predictions. 

Iqnition delay distance. The obtained water mass fraction field is given in Fig.23. It 
is seen that the flamelet model predicted self-ignition point location very close to 
the injector (x/dj-4) which is reasonably correlates with the our previous 
predictions. The inflow air temperature was high (TM500K) in conditions of Beach 
test case. So, value of N £> was high ( A^*462se C -') and changes introduced into 
the self- ignition criterion influenced results of computations very weakly for t is 
test case It is seen also (Fig.23) that maximum water mass fraction sufficiently 
rapidly trend to the near-equilibrium values downstream self-ignition point. 

The obtained Mach number field is given in Fig.24. The flowfield is slightly 
disturbed in the vicinity of the ignition point and weak compression wave is 

generated in this region. . _ llUe 

In general, reported pictures are in good agreement with our previous results. 

Reactive scalars. The comparison of the obtained reactive species mass fractions 
distributions (solid lines) with the experimental data is given in Fig.25a-d for our 
test sections where experimental data were available. The results of our previous 
predictions are presented also by dashed lines. It is seen that results of current 
predictions well correlates with both experimental data and our previous results. 
The mostly essential difference is reasonable improvement of the water peak value 
predictions in the first test section due to decreasing of the chemica 
nonequilibriumness of the detailed chemistry model compare to that used in the 

previoU Q ^ concluded that introduced modifications allowed to obtain better 

prediction of the water peak value (better that 10%) for the first test section where 
previous tests gave higher discrepancy (about 25%). The accuracy for the other test 
stations remained reasonably good. 
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TV. 3 CHENG et.al. ROUND-JET EXPERIMENT 


IV.3.1. TEST CASE 

The spontaneous Raman scattering and laser induced pre-dissociative 
fluorescence (LIPF) methods were employed by Cheng et.al. (1994) to measure 
simultaneously the temperature and concentrations of species in supersonic coaxial 
burner. The scheme of setup is sketched in Fig.26 together with inflow conditions. 
The hydrogen was injected with the Mach number M H2 =1. The injector internal 
diameter was d = 2.3610' 3 m. The vitiated air was obtained by burning of hydrogen 
in air. The high-enthalpy vitiated air was expanded through convergent-divergent 
nozzle with the Mach number M a =2 at the exit. There was small shift (*d) in 
longitudinal positioning for exit cross-sections of the hydrogen injector and the air 
nozzle. The hydrogen injector had thick lip. The lip thickness (h»0.7- 10' 3 m) was 
compatible with hydrogen jet radius. 

The flowfield was characterized by the interaction between coaxial jets with 
ambient air. So, two mixing layers were formed in the flowfield. The first one 
(where combustion mostly occurred) was between hydrogen jet and vitiated air. The 
second one was between vitiated air flow and ambient air. It was reported that the 
flame was lift-off from the nozzle lip to x/d«25 but self-ignition was asymmetric. 

The used vitiated air generator was small enough and, probably, it was 
difficult to obtain perfect premixing of contaminants after pre-burner. So, the inflow 
species and temperature distributions in the vitiated air flow were nonuniform. The 
vitiated air composition was characterized by high level of OH radicals (about 10' 3 
by volume fraction). However, independently of these disadvantages, the unique 
data about mean and fluctuating scalars in supersonic flame were obtained in this 
experiment. 

IV.3.2. COMPUTATIONS 

In the current tests, we tried to take into account presence of the ambient air 
in the computational domain. Vitiated and ambient air had different atomic 
composition. So, we applied simplified generalization of our approach to model 
presence of two mixing layers between three flows having different atomic 
composition. Two mixture fractions z\ and z 2 were introduced to describe state of 
mixing in hydrogen/vitiated air and vitiated/ambient air mixing layers respectively. 
It was expected also that combustion processes occur due to interaction of the 
hydrogen with vitiated air only. Possible role of the ambient air in chemical 
reactions was neglected. This assumption was approximately valid in the main part 
of the flame up to x/d«45 cross section as it can be concluded from the 
consideration of the data reported by Cheng et.al. (1994). 

Let us present adopted modification in more details. 

First, summary mixture fraction z was introduced in the same manner as it 
was done previously. It was defined as total mass fraction of hydrogen atoms H in 
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all species. It was varied from z = 1 in the H 2 jet to z = 0 in the ambient air, since 
ambient air did not contain hydrogen-laden species (ambient air humidity was 

neglected). The total mixture fraction had value z= z = 2 C ll2() + -^-C ()11 in 

ft 1120 Poll 

the vitiated air flow since vitiated air contained water and hydroxyl radicals with 
mass fractions C* 12() and C* )H respectively. 

The first mixture fraction z, was used to model state of mixing and 
combustion in the internal (H 2 /vitiated air) mixing layer in the same manner as it 
was done by Cheng et.al. (1994) in their processing of the experimental data. The z, 
was introduced as the normalized excess of the hydrogen atoms mass fraction z 
over its value in the vitiated air flow: 


where H is the Heaviside step function, which was introduced to reject negative 
values. The first mixture fraction z, equals 1 in the pure fuel and it equals 0 in 
both vitiated and ambient air. 

The second mixture fraction z 2 was introduced to model state of mixing 
between vitiated and ambient air in the external mixing layer. It was defined as: 

z 2 =^ + (\-~)H{z-z) 

<T T 

X, X. 

The second mixture fraction z 2 equals 1 for both hydrogen jet and vitiated air and 
it equals 0 for ambient air. 

It was expected that instantaneous values of Z\ and z 2 satisfy diffusivity 
equation in form (1.1). So, conventional Favre-averaged form of the turbulence 
modeling equations was used to calculate their mean (z,,z, ) and variance (a, 2 , a;) 
values: 

V(pVz, ) = Vp(v, / Sc, + v / Sc)Vz, ; 

V(pVa 2 ) = Vp(v, / Sc + v/Sc)Va 2 +2p^(V2,) 2 -2pN,; 

Sc, 

etc. 

These four additional equations were solved together with the main system of 
the hydrodynamics equations and turbulence model. The obtained mean and 
variance values were used for calculations of mixture fractions PDF's in internal 
and external mixing layers. The reactive scalars distributions were calculated from 
the flamelet equations which were solved vs Z\. The mixing solution vs z 2 was used 
for calculation of reactive species in the external mixing layer. The averaging 
procedure was applied in CFD for both internal and external mixing layers using 
mixture fractions PDF's. All other details of computations were not differed from 
those reported for other considered test cases. 
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The flow field was expected axisymmetric. The small longitudinal shift in 
positioning of the hydrogen injector and air nozzle exit cross sections was ignored. 
It was expected that they were positioned in the same cross section. The 
calculations were done for the same as those for the Beach test case set of the 
laminar and turbulent Prandtl and Schmidt numbers (Pr = Sc = 0.72; Pr t = Sc t =0.7; 
c CT = 1/Sc t ). 


Computational domain. The scheme of computational domain is presented in Fig. 27. 
The exit cross-section of the vitiated air nozzle was chosen as left boundary of 
computational domain. The exit computational boundary was posed in cross-section 
x/d = 43.3. The lower computational boundary was at the axis of symmetry of the 
burner. The upper computational boundary was located in ambient air. The 
transverse dimension of computational domain was equal to y/d=12 in the initial 
cross-section and it increased linearly downstream to take into account the growth 
of the external mixing layer between vitiated and ambient air. 

Boundary conditions. The left boundary included hydrogen jet, lip surface, vitiated 
air stream and ambient air. The parameters distributions in the vitiated air flow and 
in hydrogen jet were supposed to be know. The lip surface was expected adiabatic 
with no-slip velocity conditions and zero values of turbulent viscosity and turbulent 
kinetic energy. 

The flow parameters at the ambient air part of the left boundary were 
calculated using the left Riemann invariant, defined from the computational 
domain, and known total parameters for the undisturbed ambient air. The vertical 
velocity component was supposed to be equal to zero at this part of the left 
boundary. The same procedure was applied for flow parameters calculation at the 
upper boundary of computational domain 

The drift boundary conditions were used for the supersonic part of the right 
computational boundary. The condition of given pressure was used for subsonic 
part of this boundary. It was expected that pressure in this part of the exit 
boundary was 1% less than that in undisturbed ambient air. 

Inflow parameters distributions. The parameters distributions (normalized by mass- 
averaged velocity of the hydrogen jet U H2 and internal injector diameter d) which 
were used as inflow conditions at the left boundary of computational domain are 
presented in Fig. 28. They were obtained by the following manner. 

The internal flow in the convergent - divergent nozzle of the vitiated air pre- 
burner was calculated using our NS CFD code up to the nozzle exit cross-section. 
The nozzle geometry was taken from the sketch presented in Cheng et.al.( 1994). 
The total parameters after pre-heater were taken in accordance with specification 
presented in their paper also. The calculated vitiated air flow parameters in the core 
flow were very close to that reported by Cheng et.al. (1994) in their specification of 
the test case, except the vitiated air was just slightly closer to the fully expanded 
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conditions. The obtained, at the exit of the nozzle, parameters distributions were 
used as inflow profiles in the vitiated air flow. 

The fully developed turbulent pipe flow was assumed in the hydrogen 
injector. It was expected that turbulent boundary layers thickness was equal to 
radius of the tube. The longitudinal velocity, turbulent viscosity and turbulent 
energy distributions were calculated using procedure outlined in Appendix B. The 
pressure and density were supposed to be constant across hydrogen jet. The vertical 
velocity component was equal to zero in the injection tube. 

Computational grids and convergence. The overlapping grids were used to increase 
calculations accuracy near the injector lip surface. The first domain (0<x/d<3) 
included 120x140 grid nodes and the second domain (2<x/d<43.3) included 100x140 
grid nodes. The computational grids are presented in Figs. 29, 30. The boundary 
conditions at the right boundary of the first computational domain (x/d = 3) were 
formulated as it was discussed above for the right computational boundary of the 
whole computational domain. 

The CFD computations of this test case required about 3500 iterations and 15 
hours of our HP workstation. The behavior of the residual L2 norm vs iterations 
number is presented in Fig.31. Sufficiently slow convergence was due to the 
presence of very low-speed ambient air in the external part of the computational 
domain. After the first 1000 iterations, the main residual was located in these 
regions and further decreased very slowly. 

IV.3.3. RESULTS and DISCUSSION 

Flow field patterns. The obtained in computations static pressure contours and 
Mach number field are presented in Fig. 29 for the region near injector. It is seen 
series of expansion and compression waves generated in the flows of hydrogen jet 
and vitiated air around injector lip (its height was compatible with the hydrogen 
injector diameter). 

The Mach number field for the far domain is given in Fig. 32. It is clear seen 
development of two mixing layers. The inner layer (between H 2 and vitiated air) is 
reactive and the outer layer corresponds to the mixing between vitiated air and 
ambient air. It was found that the inner flow alternates periodically from being 
supersonic to subsonic before becoming supersonic again. The similar complex 
flowfield structure was reported previously by CFD modellers who studied this test 
case using different turbulent combustion models (Eklund, Drummond, Hassan, 1990; 
Hsu, Raju, Norris, 1994). Qualitative remarks regarding existence of the wave 
structure were presented by authors of the experiment also. 

Ignition delay distance. The flamelet predicted self-ignition point location in x/dM2 
cross-section. The authors of the experiment reported value x/d*25 as the 
estimation for ignition delay distance based on the long exposure visual photograph 
of the flame. From our point of view, such difference in predictions and 
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observations is explained by asymmetric ignition which was observed in the 
experiment. We expect also that the inhomogenity of composition and temperature 
of the vitiated air flow was the mostly probable cause responsible for this 
asymmetry. The inflow vitiated air inhomogenity is well-seen in Figs.33f,34 for the 
nearest test station x/d = 0.85. For example, the measured in the inflow OH mole 
fraction (Fig.33f) decreased, approximately, 3-4 times from positive (over y-axis) to 
the negative side of the vitiated air jet. We were not able to take into account all 
details of such inflow composition nonuniformity in our current axisymmetrical 
CFD. So, OH concentration (»10‘ 3 by mole fraction) was chosen in the inflow 
conditions as its typical value observed in the experiment at the positive side of the 
vitiated air jet. Thus, our predictions better reproduced behavior of the positive side 
of the jet. The temperature and H 2 0 flame peaks in experimental data are well-seen 
on the positive side of the jet even in x/d = 21.5 cross section while they are absent 
on the negative side of the jet in this cross section (please see Figs. 33d and 34). 

We performed computations where we decreased, approximately by an order, 
inflow OH concentration to its typical value for the negative side of the vitiated air 
jet. This resulted in ignition point displacement more close to its position reported 
by Cheng et.al. (1994). Such sensitivity of ignition delay distance prediction to the 
inflow radicals concentrations was mentioned for this test case by another modeliers 
also. 

Mean reactive scalars and temperature. The results of predictions for mean values 
of mixture fraction and species concentrations are given in Figs.33a-f. Our 
predictions for the positive side of the jet are re-plotted for its negative side also for 
more clearness. 

The radial distributions of the mixture fraction are presented in Fig. 33a. It is 
seen that calculations correlate fairly well with experimental data for the first three 
test sections (x/d = 0.85, 10.8 and 21.5). However, it is seen that turbulence model 
somewhat underpredicted mixing rate near the axis of symmetry for the far field of 
the jet. The same tendency was found in H 2 profiles as it is seen from Fig. 33b. We 
expect that such difference can be explained by the aforementioned asymmetry of 
the self-ignition in the experiment which could provide increasing of the 3-D effects 
in the internal mixing layer downstream self-ignition region. 

Predictions of the N 2 mole fractions were well-correlated with the 
experimental data (Fig. 33c) for all test sections. 

The obtained H 2 0 distributions are presented in Fig. 33d. We expected the 
obtained correlation as fairly good since reported CFD was done for the inflow 
conditions which approximately reproduced positive side of the jet. It is seen that, 
according to the experimental data, the positive side of the jet is within combustion 
regime, at least, from x/d<21. Only mixing was observed in experiment at the 
negative side of the jet in both 21.5 and 32.3 cross sections. As it was discussed 
previously, it was possible to reproduce such behavior of the negative side of the jet 
by decreasing of the OH inflow concentration. In any case, flamelet model 
predictions reasonably correlated with the experimental data except only slightly 
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underpredicted peak value for x/d = 32.3 cross section. However, measured for this 
test section H 2 O peak value is somewhat questionable taking into account reported 
by Cheng et.al. data on instantaneous realizations. It is seen (Fig. 38b) that H 2 O 
instantaneous data collapsed to the peak value which is about 0.6-0.65. These values 
are approximately 15% higher than that for chemically equilibrium limit. Such effect 
did not observed in another test stations and probably can be attributed to some 
instrumentation effects.. 

The reasonable correlation between experimental data and predictions was 
obtained for O 2 mole fraction in the inner mixing layer as it is reported in Fig.33e. 
The discrepancy between data and predictions in the outer mixing layer is 
explained by significant non-uniformity of the inflow 0 2 distribution provided by 
the vitiated air generator (please see data for x/d = 0.85 cross section in Fig.33e). 
This nonuniformity was ignored in current CFD. 

The calculated distributions of the OH mole fraction are presented in 
(Fig.33f). It is seen that both peak values and width of the profile were reproduced 
with the accuracy about 20-30% for the right hand side of the flame in the region 
(x/d>20). We expected this correlation as reasonable taking into account that 
reported accuracy of OH measurements was about 13%. The correlation was not so 

good upstream the ignition point (x/d<10.8). It can be explained by both 
aforementioned inflow nonuniformity and existence of slow chemical processes in 
the induction region which were neglected in our simplified criterion for self- 
ignition. 

The performed tests show reasonable correlation between measured and 
predicted by the flamelet model static temperature distributions in the combustion 
zones and internal mixing layer as it is reported in Fig. 34. Here, averaged solution 
of flamelet equations for static temperature is plotted together with experimental 
data. The obtained correlation indicates that adopted in the flamelet calculations 
simplified treatment for conditionally averaged flow velocity V = V(z) (please see 
Sec. I) operated reasonably well and it did not introduce significant inaccuracy in 
the static temperature distributions even in the case where Mach number was 
varied in wide range from 0 (ambient air) to 2 (vitiated air flow). Nevertheless, it 
should be mentioned some underprediction of the static temperature distributions 
for the external mixing layer in the final test sections (x/d = 32.3 and 43.1). We 
hope that it is not the property of our model since the same discrepancy was found 
in the Evolved and Assumed PDF modeling reported by Hsu, Raju, Norris, (1994); 
Baurle, Hsu, Hassan, (1994) and other modellers. Followed Hsu, Raju, Norris, (1994), 
it can be expected that such discrepancy could be attributed to some 
instrumentation inaccuracy for measurements in low temperature regions reported 
by Cheng et.al. ( 1994). 

Reactive scalars RMS. The obtained distributions of the mixture fraction rms are 
presented in Fig. 35 together with the experimental data. It is seen that obtained 
distributions well correlated with experimental data except only x/d=10.8 cross 
section where computations somehow overpredicted mixture fraction rms. 
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The flamelet predictions for reactive scalars rms are reported for x/d-10.8 

and 43.1 cross sections in Figs. 36 and 37. 

It is seen from consideration of Fig. 36 that rms distributions for reactive 
scalars correlate in main features with experimental data in x/d=10.8. Some 
difference between predictions and data can be explained taking into account both 
adopted in modeling simplifications and specific features of the test case. For 
example, it is seen difference in predictions of rms fluctuations in regions near 
y/ ( j = 2.5 and y/d = -2.5. These regions corresponds to the vitiated air flow. The 
nature of the experimentally observed fluctuations is the inhomogenity of the 
vitiated air composition provided by pre-heater. We did not take into account this 
feature of the inflow conditions in our modeling. Followed Hsu, Raju, Norris, (1994), 
we expect that experimentally observed values 0.01 for 0 2 and 0.04 for N 2 in the 
outer regions of the external mixing layer could be attributed to the background 
fluctuations since the data for the other test stations indicate the same level for 
these fluctuations. We expect that underprediction of the OH radical rms value was 
due to very simplified treatment of the ignition phenomenon in the flamelet 
approach. If one take into account aforementioned features, the flamelet model 
predictions for reactive scalars rms reasonably correlates in this cross section with 
results of PDF modeling reported by Hsu, Raju, Norris, (1994) except some 
underprediction of the peak rms temperature near the axis of symmetry. 

It is much more difficult to present conclusions regarding quality of the rms 
predictions for x/d = 43.1 cross section. The peak values are predicted reasonably 
well, but displacement of their locations could be attributed to the flow asymmetry 
in the experiment. 

In general, flamelet modeling provided predictions of the reactive scalars rms 
compatible with that reported for other approaches. 

Instantaneous patterns for reactive scalars. The authors of the experiment reported 
data regarding instantaneous realizations for the reactive scalars measured by their 
Raman system. The original data of Cheng et.al. are presented in Figs.38a,b by 
green points. So, we tried to use their data for simplified estimations of realisticity 
of the instantaneous reactive scalars distributions provided by the flamelet 
approach. 

It was expected that scattering of instantaneous values of reactive scalars (at 
fixed value of mixture fraction) was due to the scalar dissipation fluctuations 
(Bilger, 1989). The increasing of the scalar dissipation leads to the growth of the 
chemical nonequilibriumness in the reacting zones due to increasing of the 
reagents fluxes. Otherwise, its decreasing equilibrates the state of combustion 
chemistry. So, fluctuations of the scalar dissipation would lead to arising of different 
instantaneous patterns of the reactive scalars in the reaction zones and hence, the 
scatter plots can be used for indirect estimation of the realisticity of instantaneous 
solutions provided by the combustion model. 

Unfortunately, statistical properties of the scalar dissipation are not so well- 
investigated now as those for mixture fraction. However, both experimental data 
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(Namazicm et.al, 1988) and theoretical estimations (Monin, Yaglom, 1965) show that 
its PDF can be approximately modeled using so-called log-norm low: 

Q(A r ) = ^ exp[ -r(ln(A / N) + a 2 /2) 2 ] (IV.l) 

(2 n) N<j 2 cr" 

where a 2 ^0.51n(L t /g), L t is integral length scale and r| is Kolmogorov scale. Using 
this PDF model, it is possible to estimate upper and lower values of scalar 
dissipation which corresponds to the given probability of expectations. We choose 
the range 99% probability and calculated N max and N inm _values using log-norm 
PDF distribution for scalar dissipation. The mean value of N in relation (IV.l) was 
chosen to be equal to the value of conditionally averaged scalar dissipation at the 
stoichiometric surface provided by our CFD. The data on integral and 

Kolmogorov turbulence scales (L t and g) were taken from experimental estimations 
reported by Cheng et.al. (1994). Further, the upper and lower possible flamelet 
distributions, predicted by the model in the range of N s variation N mm SN s< CN max , 
were taken from the flamelet library. The comparison of such upper and lower 
limits for possible instantaneous solutions with the scatter data and conditionally 
averaged flamelet solution is presented in Fig. 38a, b for point y/d=l.l in x/d = 32.3 
cross section (Nf = nsec" 1 ). It is seen that upper and lower boundaries for 
temperature and main species are within scatter distributions and they are located 
close to the upper and lower boundaries of scatter plots. The conditionally averaged 
solutions are located close to the regions where scatter points clustered except H 2 0 
distribution. However, as it was discussed earlier, the observation of scatter points 
with the H 2 0 mole fraction about 0.6-0.65 is questionable. These values are 
approximately 15% higher than equilibrium chemistry peak value for H 2 0. 
Additionally, this feature was not observed both upstream and downstream 
considered test section. So, probably, it would be attributed to some 
instrumentation uncertainty. 

The flamelet somewhat underpredicted possible peak distribution of the OH 
radical as it is shown in Fig.38bT This means that model slightly underpredicted 
maximum nonequilibriumness in this point of the flame. Followed Cheng 
et.al. (1994), we can expect that this underpredictions could be due to neglecting of 
the velocity fluctuations in the flamelet calculations. The velocity fluctuations can 
cause additional temperature fluctuations and, probably, influence combustion 
chemistry in supersonic flames. 

The same kind of estimations was made for another test sections where 
scatter plots were available. They gave basically the same results. So we expect that 
these indirect estimations indicate that flamelet model correctly reproduced main 
features of the combustion processes for the supersonic flames at moderate Mach 
numbers. We expect that simplified treatment of the kinetic energy term adopted in 
the current flamelet calculations did not lead to any significant errors in 


') The maximum OH distribution corresponded to N s «17sec _l so conditionally averaged and 
maximum OH distributions were very close one to the other. 
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predictions. However, the modeling of supersonic flames at higher Mach numbers 
can require more accurate treatment for this term. 

Based on obtained results, it was concluded that flamelet approach gave 
satisfactory results for this test case. Specific features of this experiment 
(asymmetric ignition, imperfect premixing in vitiated air flow and presence of the 
ambient air) required us to introduce additional assumptions for its flamele 
modeling. Nevertheless, both temperature and reactive scalar distributions in the 
reaction zones were predicted reasonably good. The correlation for reactive scalars 
rms was somewhat worse but accuracy of the flamelet predictions was compatible 
with that reported for Assumed and Evolved PDF modeling. Additional indirect 
estimations demonstrated that instantaneous flamelet solutions were within scatter 
distributions provided by the Raman measurements. Some underprediction of the 
upper limit for OH radical instantaneous distributions can indicate possible 
necessity to use more accurate treatment of the kinetic energy term in flamelet 
equations. However, we expect that such modifications can became important for 
modeling of the supersonic flames at much higher Mach numbers than that 
considered here. 

It should be noted also, that we are not fully satisfied by the pure 
computational results of the performed simulation. Unfortunately, CFD of this test 
case were very expansive and we were practically at the limit of our HP workstation 
capabilities. We attributed this feature to the presence of the low-speed air and 
somewhat simplified formulations of the boundary conditions at the ambient air 
part of the computational domain. It is possible to hope, that the influence of these 
assumed conditions appears essential in the external part of the shear layer between 
vitiated air and ambient air and it did not influenced region of hydrogen - vitiated 
air interaction. Nevertheless, we expect to reconsider the aforementioned numerical 
aspects of this test configuration in further investigations to obtain less-expensive 
and correct modeling. 
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TV . 4 CHANG et.al. PLANAR MIXING LAYER 


IV.4.1. TEST CASE 

Scheme of setup and test conditions of Chang et al. (1993) experiment are 
given in Fig. 39. The test section was rectangular duct having the constant width 
(0.2m). Air and fuel streams converged at the splitter plate tip with a 6 degree 
convergence angle. Both air and fuel ducts were 0.05m height at the splitter tip. 
Non-reacting and reacting flow regimes were studied. Lower and upper duct walls 
were parallel in the non-reacting case. They were 1.3 degree diverged in the 
reacting case. Velocity profiles and turbulence characteristics were measured in 
several cross sections using LDV. Temperature measurements were done in cross 
sections 0.006, 0.15 and 0.3 m for the reacting case. The hydrogen-fueled torch was 
used for the flame stabilization in the reacting case. It was located in the air 
supplying duct upstream the splitter plate. Ignition in this case was observed in a 
region close to the splitter plate tip, as it can be seen from OH emission images, 
presented by Chang et al. (1995). 

VI.4.2. COMPUTATIONS 

The flow field was expected 2-D. The role of the boundary layers at the 
upper and lower duct walls was neglected. The molecular viscosity and diffusivity 
were estimated based on air molecular diffusivity and fixed laminar Prandtl and 
Schmidt numbers Pr = Sc = 0.72. The computations were done for values of Sc, 
= l/c CT = 0.5. 

Probability of z = 1 value was non-negligible along the whole flowfield for the 
considered mixing layer case. So, small modification of the PDF model was applied 
to take into account flow intermittency effects at fuel-rich edge of the mixing layer 
(z= 1). Followed Kuznetsov and Sabelnikov (1990), the mixture fraction PDF was 
expected to be symmetrical with respect to z = 0.5. Probability to observe both z = 0 
and z = 1 values in the same point of the mixing layer was neglected. Based on such 
assumptions, mixture fraction PDF was modified to the form suggested in Kuznetsov 
and Sabelnikov, (1990). In regions where z < 0.5 , conventional relations of Sec. I were 
applied for calculations of mixture fraction PDF and intermittency factor: 

./(z) = (l -y)6(z) + yP l (z); y = min|l.31z 2 / (z 2 +cT ); lj . 

In regions where z > 0.5 , the mixture fraction PDF and intermittency factor were 
calculated as: 

;y(z) = (1 - y, )8(1 - z) + y ,P t (l - z) ; y, = min{l 3 1(1 - z) 2 /((1-z) 2 +cr); l} . 

The flamelet model calculations showed that was very small (N ( J <10 3 

sec 1 ) for conditions of this test case. It indicated that self-ignition was impossible 
inside the test section and it was necessary to apply external device for 
initialization of combustion wave. This qualitative prediction correlates with the 
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necessity to use hydrogen-fueled torch in the experiment for flame stabilization 
inside test section. So, it was expected that flame was stabilized due to the presence 
of the torch. That is why, flamelet library solutions up to the extinction critical 
value of scalar dissipation N ( c ' } were used. All other features connected with the 
presence of torch (some temperature rise in the part of air stream; presence of small 
amount of combustion products) were not taken into account. 

The scheme of computational domain is given in Fig. 40. The left boundary 
( X = 0) was located near the splitter tip cross-section. Right boundary was set at 
x = 0.4m, which is 0.07m downstream the last presented test section. Conditions at 
the inlet and exit boundaries were posed in the same way as it is often done for 
subsonic flows. Velocity and static temperature profiles, profiles of turbulent 
characteristics and mixture fraction were supposed to be known at the inlet 
boundary while the pressure was extrapolated to the inlet boundary from the 
computational domain. Fixed pressure condition (0.1 MPa) was applied at the exit 
boundary. Normal derivatives of all other parameters at the exit plane were 
determined from computational domain. The upper and lower walls were 
considered as inviscid with zero normal velocity component. All normal derivatives, 
which were needed to estimate viscous stresses and corresponding diffusion fluxes 
on the wall, were equal to zero. 

The parameters distributions at the inlet boundary were taken from LDV data 
at the first test section reported by Chang et.al. (1993). Streamwise velocity in core 
flows was set equal to that measured in the experiment, i.e. 394 m/ s at the air side 
and 134 m/s at the fuel side for non-reacting case; 390 m/s and 137 m/s 
respectively for reacting case. Static temperature in core flows was estimated using 
measured values of total temperature and Mach number. Approximate analytical fit 
to the experimental data was used to model the inlet profile of the cross-stream 
velocity component. The thickness of the boundary layers on the lower and upper 
surfaces of the splitter plate was adjusted to the measured streamwise velocity 
profile. The eddy viscosity and streamwise velocity profiles inside the mixing layer 
were obtained using the procedure described in Appendix B and they were slightly 
smoothed near splitter lip. The value of turbulent kinetic energy in the core flows 
was estimated using experimental data on velocity components rms as K«u’ 2 /2 + v' 2 . 
In the non-reacting case, inlet turbulent kinetic energy profile in the mixing layer 
was obtained using procedure of Appendix B. In the reacting case, turbulent kinetic 
energy was approximated using experimental data on velocity components rms as 
K~u' 2 /2 + v' 2 . The uniform step initial profile was expected for mixture fraction z. 
The mixture fraction variance o 2 was expected to be zero at the entry boundary. 

To estimate the value of turbulent viscosity in the core flows, an approximate 
relation, which is often used in turbulence modeling, was used: Vi^O.lu'L, where L, 
is the turbulence integral length scale. Unfortunately, no data were available on 
turbulence scale in this flow, so the value of L t was approximately estimated as 0.1 
of the supplying duct height, i.e. L^~5 mm. As it was found in computations (see 


39 



below), the flowfield was sensitive to the inlet value of L t , especially in the reacting 

The non-dimensional initial distributions of the streamwise and cross-stream 
components of the flow velocity U = — and V = — ; eddy viscosity v t - — and 


turbulent kinetic energy K = are given in Fig.41 where u a is the air flow velocity 

K 

in the core flow and h is the test section height at the splitter tip cross-section 
(h = 0.1m). 

Calculations were done using the nonuniform grid with 100 cells in 
longitudinal direction and 80 cells in transversal direction (see Fig.42). Grid was 
clustered to the duct centerline and to the inlet plane. 

The convergence was estimated by the L 2 norm for the residual of continuity 
equation. The L 2 norm decreased 4 order during 800 iterations which required 
about 5 hours of Pentium-200 PC. 


VI.4.3. RESULTS and DISCUSSION 


Non-reacting case. The obtained profiles of mean streamwise velocity are presented 
in Fig. 43 for eight cross-sections where experimental data were available. It can be 
seen that streamwise velocity profiles were in good agreement with the 
experimental data. The mean cross-stream velocity distributions were not predicted 
so accurately (see Fig.44). The calculated cross-stream velocity distribution rapidly 
decayed downstream the splitter, while experimental data showed existence of its 
sufficiently high negative values far downstream the splitter plate. 

Profiles of calculated turbulence kinetic energy are compared to that, 
obtained from experimental data in Fig. 45. The relation K»u' 2 /2 + v 2 was used here 
to estimate turbulence kinetic energy from experimental data on rms velocity 
components. The agreement is reasonably good for the first four test stations 
(x< 150mm). Some disagreement in turbulent kinetic energy distributions was found 
far downstream splitter plate (x = 300 and 330mm test sections). However, as a 
whole, predictions of the turbulence kinetic energy distributions were acceptable 
taking into account simplified relations for its estimation from experimental data. 

The significant disagreement between calculated and measured Reynolds 
stress profiles was obtained. The Reynolds stress provided by our CFD was obtained 


using relation (u'v') = - 


du dv^) 

h 

dy dxj 


. The computations 


overpredicted experimental 


data several times (usually 5-7) as it is reported by Fig.46 where the calculated and 
measured Reynolds stress normalized by the slip speed square vs. self-similar 
coordinate q = y/x are shown for x = 300mm test station. The same result was found 
for other test sections also. This result was unexpected since we sufficiently 
accurately predicted distributions of the mean streamwise velocity component. So 


40 



we compared results of our CFD with data regarding (u'v') behavior in planar 
mixing layers obtained in different experiments and summarized in Abramovich 
et.al.f 1984). It is seen (Fig. 46) that our results correlate with these data reasonably 

well. 

The obtained discrepancy can not be explained by the influence of the 
discretization as it is reported by Figs.47a,b where results obtained for three 
different grids are shown (100x80, 150x120 and 200x160). Followed Lai and Raja 
(1993), who performed CFD modeling of this experiment, we can conclude that 
even results presented for "rough'' 100x80 grid mesh provided good resolution for 
the parameters distributions in the mixing layer. 

We expect that the obtained result indicates that even non-reacting 
counterpart of the Chang et.al.(1993) experiment contained some abnormal features 
of the turbulence spectra which could not be reproduced using conventional 
turbulence modeling relations. 

Reacting case. The difference between predictions and experimental data was found 
for the reacting counterpart of the Chang et.al. experiment also. The calculated 
streamwise velocity distributions are compared with the experimental data in Fig. 48. 
It is seen that our results correlate with the experimental data in region x< 150mm. 
Significant difference between results of computations and data was found for the 
far field (x = 300 and 330mm cross-sections). The experiment indicates a much 
significant penetration of the mixing layer into the upper fuel-laden stream. It is 
seen the underprediction of the mixing layer growth obtained in computations. 
Again, as in the non-reacting case, it was impossible to reproduce correctly cross- 
steam velocity distributions as it is reported in Fig. 49. The maximum in the 
turbulent kinetic energy distributions was predicted more correctly but significant 
underprediction of the mixing layer width and its positioning was found in the far 
field also (Fig. 50). 

As for the non-reaction case, it was found that grid influence could not 
explain obtained discrepancy as it is reported by Fig. 51a, b. It is seen that the use of 
finer grid influenced the results of computations very slightly. 

One could expect that operation of the torch ignitor in the air supplied duct 
could lead to the development of the large scale fluctuations in the reacting case. 
So, we tried to estimate possible influence of such large-scale fluctuations on the 
results of predictions. For this purpose, the turbulent viscosity in the incoming air 
flow was increased ten times which corresponded to the increasing of the integral 
length scale approximately up to the value compatible with the air supply duct 
height (L t «50mm). Such modification provided much better correlation of both 
width and location of the mixing layer in the far field as it is reported in Fig.52. 
The cross-stream velocity distributions became slightly closer to the experimental 
data also (Fig. 53). However, both abnormally high value of the integral length scale 
which was necessary to obtain such improvement and absence of data about its 
behavior in the experiment, does not allow to consider obtained result as any kind 
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of final predictions. It should be considered as rough estimation which indirectly 
indicates possible role of some unconventional scales distributions in the turbulence 
spectra or presence of some large-scale fluctuations in the experiment. 

The predicted by the model mean temperature distributions are presented in 
Pig 54 for x = 300mm cross section. It is seen that the maximum in the temperature 
distribution is reproduced by the model but pure accuracy in flowfield predictions 
lead to pure accuracy in predictions of both width of the temperature distribution 
and positioning of the mixing layer. The estimation provided by the large-L t 
computations is presented also. It is seen that such modification somewhat improve 
predictions but, as it was discussed previously, it can not be considered as the final 
result since the presence of the turbulent fluctuations having the integral length 
scale compatible with the duct height is quite questionable. 

As a whole, results of the current calculations correlates with those reported 
by Lai and Raja (1993), who studied capabilities of their RPLUS and composite PDF 
models for this test case. As the aforementioned modeliers, we obtained the same 
discrepancy in predictions. We agree with them that presence of the hydrogen 
torch provided some element of uncertainty in inflow conditions for the considered 
test case. Such feature could be responsible, in particular, for discrepancy between 
results of predictions and experimental data. However, it should be mentioned that 
we found difference for Reynolds stress distributions even in non-reacting case. This 
finding makes this test case different from typical ones. We hope that obtained 
discrepancy can not be attributed directly to the flamelet model since: 

(i) discrepancy was observed even for nonreacting counterpart of the experiment; 

(ii) the same discrepancy was obtained by Lai and Raju (1993), who studied this test 
case using other turbulent combustion models. 

We expect that performed estimations showed some abnormal behavior of the 
turbulence spectra or presence of some large-scale fluctuations in the experimental 
setup. The observations of the authors of the experiment indicated presence of the 
"large-scale shear layer flapping motion" in the test section. Probably, better 
correlation could require to model more accurately features of inflow and outflow 
ducts of this experimental setup. Unfortunately, in spite of the good contact with 
one of the experiment authors (Dr.C. Chang, LeRC), causes of such behavior of the 
mixing layer are not quite clear for us yet. This require us to retain further 
consideration of this test case for the next investigations. 
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V. MODELING TESTS OF MODIFIED FLAMELET APPROACH FOR FLAMES 
WITH LARGE IGNITION DELAY DISTANCE 

V.l. MODIFIED FLAMELET EQUATIONS (MFL) 

In the previous tests, including those presented in Sec. IV of the current 
Report, we used the following form of the flamelet equations: 


N f 


d 2 C 


dz 2 


— + R„ =0 


a= 


(V.l) 


where C u is mass fraction of a-specie; z is mixture fraction; N s is scalar dissipation 
at the flame front; R„ is chemical production term. All quantities in (V.l) are 
considered as conditionally averaged over mixture fraction z. Such form of flamelet 
equations is obtained from the instantaneous conservation equations for reactive 
scalars by neglecting the convective terms and using the mixture fraction as an 
independent variable instead of space coordinate. The account, presented by Peters 
(1984), Bilger (1982) and Kuznetsov (1982), showed that, formally, such equations 
are valid in the regions where thickness of the reaction zones is small enough. As a 
rule, such situation is observed in flames far downstream ignition point. These 
equations are approximately valid in the mixing region located upstream ignition 
point due to similarity of the scalars transport equations when the difference in 
molecular diffusivities is ignored. However, the rejected terms can play important 
role in the ignition region where transition from mixing to diffusion flame 
combustion regime is occurred and reaction zone can be sufficiently thick. 

The improper behavior of the flamelet model in the vicinity of self- ignition 
point significantly restricts its capabilities for lift-off diffusion flames. It is a serious 
disadvantage for the modeling of the high-enthalpy compressible flames since, in 
many cases of such flames, test geometry does not contain flameholders for reliable 
flame stabilization near injector and flame self-ignite far downstream injection 
point. 

Trying to improve flamelet behavior in self-ignition region, we introduce into 
our flamelet equations convective term. For the first approximation, the longitudinal 
convection was expected to be dominant, since analysis presented by Li, Bilger 
(1993) and Klimenko (1995) showed that variation of conditionally averaged values 
across mixing layer is not too significant as a rule. As the result, modified flamelet 
equations had the form close to the simplified form of the Conditional Moment 
Closure (CMC) equations proposed by Bilger (1993) and Klimenko) 1990): 


V. 


DC. 


N 


d 2 c 


-+R 


(V.2) 


dx dz " 

where V z is flow velocity conditionally averaged over z; x is space coordinate. 

It is necessary to note that, even in the steady-state case, full CMC equations 
are much more complex and 4-D (3 space coordinates + mixture fraction z are used 
as independent variables). Their accurate realization in CFD is too computationally 
expensive procedure. So, the goal of our current modeling tests was not to test 
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capabilities of CMC modeling. Our studies were confined by searching of minor 
modifications of the flamelet model which can provide its more reliable operation in 
near-ignition regions without much complication of the approach as a whole. That 
is why, the form of eqs.(V.2) was chosen as possible candidate. 

The simplest closure hypothesis were applied in the current tests. It was 
expected that the conditionally averaged flow velocity V z equals to the mean value 
of the streamwise velocity component at the stoichiometric surface. 

V, - u | n 

The scalar dissipation N and mixture fraction z were expected statistically 
independent in fully turbulent regions. Conditionally averaged value of scalar 
dissipation at the stoichiometric surface was used for approximation of N in 
eqs.(V.2) in the same manner as it was done in FL approach (please, see eq.1.8): 

nI 


^ ( y I are the mean values of the scalar dissipation IN and 

Z- ZS I 7. ZS 

intermittency factor y calculated under the condition that mean value of mixture 
fraction z = z s . 

In such treatment only mean flow velocity distributions and turbulence 
mixing characteristics were necessary for the modeling tests of the modified 
flamelet equations (MFL). For modeling tests of MFL approach, we extracted these 
distributions from results of our flamelet CFD reported in Sec. IV. So, inverse 
influence of the combustion model on the mean and fluctuating flow fields was 
neglected in current modeling tests. 

V.2. SOLUTION PROCEDURE 


To solve MFL eqs. (V.2), we slightly modified our flamelet solver FLSLV described 
in Bezgin et.al. (1996). The approximation of the derivatives remained the same, i.e. 
second order central difference for the diffusive term and first order upwind 
difference for convective term. The treatment of chemical source term was not 


changed also, i.e Jacobian 


dK r 


DC 


pz 


was used to provide implicit treatment of the 


source term. The resulting algebraic system was solved using matrix sweeping at 
every level in marching direction x. The accuracy of the solution was provided by 
the small limiting value of concentrations increment at each marching step. Usually 
the increment of all concentrations was set to be less than 1%. This condition was 
used for the choice of the step in marching direction. The accuracy of the solution 
was checked by the variation of the mentioned limit of the increment. 
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V.3. COMPARISON OF THE FL AND MFL SOLUTIONS BEHAVIOR. 


Conditions of the previously considered Cheng et.al.( 1994) and Burrows— 
Kurkov (1973) test cases were used for comparative modeling tests. The u|, and 

Nf distributions were taken from the results of our FL CFD reported in Sec. IV of 
the Report. 

The comparison of the FL and MFL instantaneous solutions is presented in 
Figs.55a-d where reactive scalars (H 2 O f 0 2 and O) and static temperature are 
plotted vs mixture fraction z and scalar dissipation N. Here, decreasing of the scalar 
dissipation corresponds to the increasing of the distance downstream inlet 
boundary. It is seen significant difference in properties of FL and MFL 
instantaneous solutions near self-ignition point. FL generated "instantaneous" (and 
non-physical) consumption of the partially premixed reactants which initially 
penetrated into fuel-rich and fuel-lean parts of the mixing layer upstream ignition 
point. The FL eguations predicted instantaneous formation of thin diffusion flame 
reaction zone near stoichiometric surface which is accompanied by mixing in fuel- 
rich and fuel-lean parts of the mixing layer. The MFL demonstrated more gentle 
consumption of the partially premixed reactants with relatively slow formation of 
the diffusion flame reaction zone in the vicinity of the stoichiometric surface. The 
development of the diffusion flame front was accompanied in MFL calculations by 
decaying chemical reactions between partially premixed reactants in fuel-rich and 
fuel-lean parts of the mixing layer. It is seen that behavior of the MFL and FL 
solutions became closer one to the other in the "well-developed" combustion 
regions far downstream ignition point. However, even in these regions there is some 
difference in predictions. The MFL predicted that some amount of the oxygen 
initially penetrated deeply into the fuel-rich part of the mixing layer does not react 
since temperature in these regions is small enough. Such unburned oxygen "bump" 
in fuel rich zones is well seen far downstream self-ignition point as it is reported by 
Fig. 55a. 

In general, MFL solutions demonstrated more realistic properties in self- 
ignition regions. They predicted gentle consumption of the partially premixed 
reactants and finite rate development of the diffusion flame front. Qualitatively, 
MFL predictions correlated with known observations of the diffusion flame 
formation in lifted flames with significant ignition delay distance (Kioni, et.al. 1993). 
However, MFL equations demonstrated more non-equilibrium properties of 
instantaneous solutions compare to the FL predictions downstream self-ignition 
region. The predictions of FL and MFL equations became the same only in the 
near-equilibrium region very far downstream ignition point. 

Further estimations were done for the averaged solutions of the MFL 
equations. Using mean and variance mixture fraction distributions provided by FL 
CFD, we estimated possible ignition delay distance and mean reactive scalars 
distributions for MFL equations in conditions of Burrows-Kurkov test case. It was 
found that MFL significantly overestimated ignition delay distance compare to both 


45 



experimental data and FL predictions. Only increasing of OH concentration (by an 
order from its equilibrium value which was used in FL predictions) in vitiated air 
allowed to obtain self-ignition region, in MFL estimations, approximately in the 
same location where it was predicted by the FL and observed in the experiment. So, 
quantitative correlation of the current MFL predictions with experimental data was 
sufficiently poor. But qualitative results of its predictions were reasonably realistic. 
The most interesting of them are given in Figs. 56, 57. 

The mean oxygen mass fraction behavior in the conditions of Burrows- 
Kurkov test case is illustrated by Fig. 56 where both predictions of FL and 
estimations for MFL approach are presented. It is seen that averaged MFL solution 
predicted gentle transition from mixing to burning regimes with very slow reacting 
of the oxygen initially penetrated into the fuel-rich parts of the mixing layer 
upstream ignition region. It is important to note that such oxygen “bump" exists in 
fuel-rich regions up to the exit plane of the test section. This MFL prediction 
realistically reproduces mean O 2 behavior in conditions of Burrows-Kurkov 
experiment as it is shown in Fig. 57. Here, both FL and MFL predictions for mean 
0 2 profile are given together with experimental data. It is seen the presence of 
unburned oxygen in the fuel-rich zone of the flame (y<0.01m). Thus, MFL indicated 
that elements of partially-premixed combustion could be important up to the exit 
plane of the test section. This prediction correlates with our conclusions regarding 
Burrows-Kurkov test case which were presented in Sec. IV. 

Unfortunately, quantitative predictions of MFL were sufficiently poor as it is 
demonstrated by Fig. 57 also. MFL significantly overestimated chemical 

nonequilibriumness and hence underestimated heat release in the flame. So, both 
peak value of the mean H 2 0 concentration and displacement of the mixing layer 
from the lower wall were in poor agreement with the data. 

Nevertheless, we hope that current studies toward MFL development gave 
encouraging results. Of course, quantitative predictions of MFL were sufficiently 
poor but its qualitative behavior was quite realistic. The latter result correlates with 
the data of Mell et. al. (1994) who compared FL and CMC properties with the 
results of DNS for the case of diffusion flame development in homogeneous 
decaying turbulence for the system with one-step modeling reaction. 

So, we hope that it is possible to develop non-expensive procedure based on 
minor modifications of the FL approach for much more realistic modeling of high- 
enthalpy delayed flames. We expect that the obtained discrepancy was due to the 
very simplified closure models for the scalar dissipation and conditionally averaged 
velocity which were used in current modeling tests. We hope to obtain better 
accuracy of predictions in next investigations using more justified treatments for 
these terms proposed in different studies (Kuznetsov, Sabelnikov 1990; Li,Bilger, 1994 ; 
Bilger, Klimenko, 1993). 
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SUMMARY and SUGGESTIONS FOR NEXT STUDIES 


The current CFD tests for compressible flames were done using flamelet 
approach reported previously in Bezgin et.al. ( 1996). 

To improve previously obtained inaccuracy of model predictions, the 
following modifications were introduced: 

(i) the modeling coefficients in mixture fraction variance equation were slightly 

changed within the range reported in the literature to decrease the 
mixture fraction fluctuations intensity provided by our turbulence model; 

(ii) previously neglected, selective efficiencies for the molecules acting as a third 

bodies in three-molecular reactions were introduced into the detailed 
chemistry model; 

(iii) the second critical value of the scalar dissipation N^ 1 was used in self-ignition 

criterion instead of the critical value at quench limit which was used in 
our previous predictions; 

(iv) presence of free radicals was expected in the vitiated air provided by the pre- 

burners. Their mass fractions were taken according to the experimental data 
(if available) or were estimated from equilibrium composition computations. 

All other details of current flamelet modeling were not differed from those 
reported in Bezgin et.al.f 1996). 

We considered 4 test cases which covered 3 main types of the flow field 
(round jet, planar wall jet and planar mixing layer) and two classes of flames (with 
and without large ignition delay distance). All current CFD were done using full 
system of averaged 2D Navier- Stokes equations. The Secundov's one-equation 
model "v t -90" was used for turbulence modeling. The hydrogen combustion 
chemistry was approximated by the detailed H 2 /air kinetics scheme proposed by 
C.Jachimowski (1988) which included 33 reactions between 13 species. 

ft was found that introduced modifications allowed to improve results for the 
previously considered test cases (Beach round jet and Burrows-Kurkov wall jet 
experiments) and to obtain reasonable agreement with Raman data obtained by 
Cheng et.al. for round jet supersonic flame. As a rule, peak values of H 2 0 and 
static temperature were reproduced with the accuracy within 10-15% in these test 
cases. Reasonable correlation was found for reactive scalars rms for conditions of 
Cheng et.al. test case. 

It was possible also to correlate location of the self-ignition point using 
modified self-ignition criterion with the accuracy about 20-30%. This result is much 
better than that obtained in our previous investigations (1.5-2 times). However it 
should be noted, that significant influence of the OH concentration value in the 
vitiated air flow on the prediction of self-ignition point location was found in 
current tests. This feature prohibited to obtain more accurate estimations for 
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capabilities of modified self-ignition criterion since only one of considered test 
cases (Cheng et.al. round jet experiment) contained information about OH radical 
value in vitiated air flow. However, even in this case, it was difficult to obtain more 
precise estimation due to inhomogenity of the vitiated air flow composition and 
temperature which led to asymmetric ignition in the experiment. 

In general, introduced modifications allowed to obtain predictions of the 
reactive scalars and temperature distributions compatible (sometimes even better) 
with that reported for Evolved and Assumed PDF modeling. The model 
demonstrated good results for supersonic flames with small or moderate ignition 
delay distance, even in the case when very simplified modeling relations were used 
for its closure. 

We met difficulties with modeling of Chong et.al. (1993) planar mixing layer 
experiment. But we hope that our explanations presented in Sec.IV.4 of the Report 
showed that the inaccuracy of predictions did not connected with flamelet model. 
Probably, more accurate treatment of the large-scales behavior in the turbulence 
model and/or inflow and outflow conditions of the experimental setup could 
improve accuracy of CFD predictions also. We hope to study possible directions for 
improvement in next investigations. 

However, independently on relative success of the approach, tests showed 
that current flamelet version does not provide reliable predictions of all flame 
properties in situations when ignition delay distance is very large (Burrows-Kurkov 
test case). The situations, where diffusion flame is significantly lift-off from the fuel 
injector, can occur in supersonic flame sufficiently frequently. So, improvement of 
the flamelet behavior in self-ignition regions is the key problem for its further 
generalization for supersonic flames. 

In the current study, we started development of modified flamelet approach 
(MFL) for more proper model operation in flames with large ignition delay distance. 
The convective term was introduced into the flamelet equations by analogy with the 
Conditional Moment Closure approach of Bilger-Klimenko. Modeling tests showed 
that, qualitatively, MFL predictions correlated with known observations of the 
diffusion flame formation in lifted flames with significant ignition delay distance. 
However, we obtained that MFL equations somewhat overpredicted nonequilibrium 
chemistry effects far downstream self-ignition point. So, quantitative correlation of 
the current MFL predictions with experimental data was poor. We expect that this 
was due to very simplified closure model for the conditionally averaged velocity 
which was used in current modeling tests. We belive that the obtained 
disadvantages can be improved in future by using more justified closure relations 
and it is possible to develop non-expensive procedure based on minor modifications 
of the FL approach for much more realistic modeling of high-enthalpy delayed 
flames. 
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Suggestions for next studies: 

Based on results of performed study and preliminary discussions with our supervisor 
Dr.L.Povinelli (NASA LeRC) we propose the following main directions for future 
investigations: 

1. Comparative tests of the flamelet and Evolved PDF modeling in the same 
“surrounding" conditions. 

We expect to prepare Monte-Carlo simulator for Evolved PDF modeling and to 
perform comparative tests of the flamelet and Evolved PDF for the 3-4 cases of sub- 
and supersonic H 2 /air diffusion flames. Other elements of the modeling (turbulence 
model, initial conditions, kinetics model) will be the same for both approaches in 
these tests. The results of such CFD will be used for more detailed analysis of 
similarity and difference of the approaches (joint PDF form, micromixing model, 
different high-order correlations, etc.). We expect that such comparative analysis 
would be helpful for further development of both approaches. 

2. Further improvement of the flamelet modeling. 

The current investigations showed that the incorporation of the convective term into 
the flamelet equations improves model behavior in the vicinity of the ignition point. 
However, its quantitative predictions are not satisfactory now. So, we propose to 
continue studies toward model improvement. For this purpose, we would like to test 
capabilities of different closure models for scalar dissipation and conditionally 
averaged velocity approximation. We would like also to develop approach for 
incorporation of such MFL model into CFD codes and to perform direct CFD tests. 

3. Turbulence modeling 

The accuracy of the turbulence model crucially influences results of turbulent 
combustion modeling. The main part of the current CFD tests was done using 
sufficiently old one-parametric "v t -90" turbulence model. New versions of one- and 
two-equations turbulence models have been developed recently in ECOLEN&CIAM 
by Prof. A.Secundov. They include more accurate treatment of different effects 
which are important in modeling of the compressible flames (dilatation and 
compressibility effects, difference between turbulence in axisymmetrical and planar 
flows, role of the external large-scale turbulence, etc.). We propose to study their 
possible benefits in flamelet modeling of compressible jet flames. 

We propose also to study possible benefits of the second order closure 
modeling in turbulence model equations for passive scalar since our current 
modeling was based on the simplest gradient diffusion model. 
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Appendix A. DESCRIPTION OF THE MATHEMATICAL MODEL 


The applied mathematical model consisted of the following main blocks: 

(i) hydrodynamics conservation equations; 

(ii) turbulence model; 

(iii) flamelet model; 

(iv) averaging procedure; 

(v) CFD/flamelet coupling procedure 

(vi) thermodynamics and detailed chemistry models of Appendix C. 


A.1 Hydrodynamics Model 

The Favre-averaged, steady-state conservation equations for 2-D or 
axisymmetric turbulent flows were written as: 


dF dG . R n 

— + — + j— = 0 

dx dy y 


(A.1) 


where: x,y are axes of Cartesian (j = 0) or cylindrical (j = l) coordinates. The 


fluxes and source term in (A.1) were given as: 


r pu ^ 


^ pv ^ 

pu 2 + p-T xx 


puv - X yx 


r G = 



puv - x xy 


PV + P-'tyv 

IpuH - ux xx - vx xv + q X J 


IpvH - ux yx - vx w + q yy i 


R 


pv 

puv -T yx 

PV 2 "Xyy +T 00 

pvH - dx - vx + q 


y/ 


(A.2) 


where: p is mean density; u,v are the Favre-averaged components of the 

velocity vector U ; p is pressure; H is total enthalpy; q x ,q y are heat flux 

components. The viscous stresses and heat flux components in (A.2) were 
written as: 

t„=P(v.+vX2— --V-0), t„=P(v,+vX2 — --V-0), 
dx 3 dy 3 
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(A.3) 


Txy = X yx =P(V, + V)| 



= P(v t + v)(2 — V ■ U) , 


q x = -P(v t / Pr t + v/ Pr) — , q = -p(v t / Pr t + v/ Pr) — 

ox d y 

where is turbulent viscosity and v is laminar kinematic viscosity; V • 0 is 


velocity divergence; h is static enthalpy; Pr and Pr t are laminar and 
turbulent Prandtl numbers respectively. 


A.2 Turbulence model 

The turbulent viscosity v t in (A.3) was calculated using "v r 90" 
Secundov's turbulence model: 

d 


dpuv t dpvv t + . pvv t _ d 
dx dy y dx 


P( c i v t +V )-^T 


+ - 


dx J dy 


-t \^ v t 

p( Cl v t +v — 
dy 
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p(c,v t +v) 


dv, 

dy 


(A. 4) 


+ c 2 pv t G + c 3 v t 


^_dp _dp^| 
u — + v — 
V dx dy J 


_ 2 G 2 _c s v,-+c 6 v t v 

C 4 pv t — -p- — 

a S 


where: a is speed of sound; S is minimal distance to the wall; 


c, = c. 


y 4 G : 
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V (3 0 v t ) 2 + y 4 G 2 J 


v, 2 + 1 1.2v,v + 12.8v 2 
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The Reynolds-averaged velocity components u,v which met in the 
RHS of (A.4) were modeled using gradient diffusion relations: 

u = u -i — — — ; v = v + — — — , where Sc. is turbulent Schmidt number. 

p • Sc t dx p • Sc t dy 

The turbulence kinetic energy K, mean mixture fraction z = pz / p and 
its variance cr=pz"z"/p were calculated using semi-empirical transport 
equations of the turbulence modeling: 
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dpuK dpvK . pvK _ 8 
dx dy * y dx 


p(l 4v, + v) 


dK 
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§H§) 
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- 2pN 


where mean scalar dissipation N was approximated as N = 0.1 


Ko 


A. 3 Flamelet Model 

The flamelet equations were taken in form (1.5) of Sec. I: 
d 2 C 

N s ^ + R a =0 a= 1,...,J (A.5) 

dz" 

U 2 

h = ( H h - H A )z + H A 

2 

where C u is mass fraction of a-specie; R< t is the chemical production term; J 

j 

is total number of reactive species; H is total enthalpy; h = Vh a C a is static 

a I 
T 

enthalpy; species specific enthalpies h a =JCp a dT+Ah a (T„) are used taking 

To 

into account species heats of formation Ah a at reference temperature 

T 0 =298. 15K. The parameter N s is the value of instantaneous scalar 
dissipation at the stoichiometric surface; superscripts F and A denote total 
enthalpies of fuel and oxidizer respectively. 

The kinetic energy term in eqs.(A.5) was approximated using relation: 
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U-U A .. r n 

U F - U A 

where U A ,U F are the mean flow velocities in the air and fuel core flows 
respectively; P is exponent which was chosen as P«l/Sc t (Sc t is the turbulent 
Schmidt number). 

The boundary conditions for eqs. (A.5) were posed in accordance with 
the test cases specifications at z = 1 (pure fuel) and z = 0 (pure oxidizer): 

z = 0 H = H A ;C a = C A a= (A.6) 

z=l H = H f ;C„ = c: 

where superscripts F and A denote composition and total enthalpies for the 
flows of fuel and oxidizer respectively. 

The flamelet eqs. (A.5) with boundary conditions (A.6) were solved 
parametrically over scalar dissipation N s e[0, N ( c ] r } ) and pressure in reaction 
zone p s . Here, is critical value of scalar dissipation at quench limit. The 
obtained solutions were collected into the flamelet library in parametric form 
on z, N s and pressure p s : 

C„ = C„(z, N s , p s ) a= 1,...,J (A.7) 

T = T(z, N s , p s ) 

The flamelet libraries for conditions of Beach test case (Sec.IV.2), 
Cheng et.al. test case (Sec.IV.3) and Chang et.al. test case (Sec.IV.4) were 
generated for fixed value of pressure in the reaction zone (p s = 0.1MPa). The 
flamelet library for conditions of the Burrows-Kurkov test case was generated 
in the range of pressure variation 0.08-0. 12MPa. Computational procedure 
for flamelet library generation was given in our previous Report (Bezgin 
et.al, 1996). 

Averaging of the flamelet solutions and their coupling with CFD code 
was done in accordance with items A.4 and A.5. 


A.4. Averaging Procedure 

In the current tests, averaging procedure remained the same as that 
used by us previously ( Kuznetsov , Sabelnikov, 1990). The Favre joint PDF of 
mixture fraction z and scalar dissipation i7(z,N s ) was considered in a 
following form: 

;/(z,N s ) = ^ P(z, N s ) = (1 - y)5(z)6(N s ) + yP t (z)5(N s - Nf ) 

P 
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where y is the intermittency factor; P t is the mixture fraction probability 
density function in a turbulent mixing layer; 5 is the Dirac function. 

The intermittency factor y was calculated using approximate relation: 
fl.31/(l + o 2 / z 2 ) if a/z>0.555 
1 if a/ z< 0.555 


y = 


(A.8) 


where z = pz/p is Favre averaged mixture fraction and a 2 =pz"z"/p is 
mixture fraction variance. Self-similar solutions of PDF equation were used 
to approximate the mixture fraction PDF Pt(z) in the turbulent mixing layer 
(0<z< 1): 


1404 z 

— ^ — Ai(1.788— — 2.338) jf y < 1 
z, z, r 


P,(z) 


1 




■exp 


KC 


( Z ~ Z ) 2 

2a 2 


A 


if y = 1 


(intermittent regions) 
(non - intermittent regions) 


where z t = z/y is the mixture fraction value conditionally averaged over the 
moments when the turbulent mixing layer is observed in a given point; Ai is 
Airy function. 

The conditionally averaged value of scalar dissipation at the 
stoichiometric surface was approximated as: 

- n| 

Nf =— p 5 - 2 - (A.9) 

y\-t » 

where n| , y|. ^ are the mean values of the scalar dissipation N and 


intermittency factor y calculated under the condition that mean value of 
mixture fraction z = z s . 


A.5. Coupling with CFD solver. 


The following procedure was used for the flamelet library 
incorporation into the compressible flow hydrodynamics solver. The 
"effective" heat capacities CP a of the reactive species were introduced in the 
same manner as it was done by Gaffney et.al. (1992): 


CP„ = jc„dT/(T-T„) 


(A. 10) 


where T 0 is the reference temperature. 

Using the thermal equation of state for mixture P = pRT/p and eqs.(A. 10) one 
can obtain the following "effective" form for the total enthalpy H: 


H 


(r-i)p 


p + Q + <^ 


(A. 11) 
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where: 


1 



R 

CP • [i 


is ‘'effective" heat capacities ratio; 


Q = (^C a Ah (1 - CP T 0 ) is "effective" heat of mixture formation; 


g is the mixture molar weight; Ah„ are species heats of formation and 
R = 8.31 J/(mol K) is universal gas constant. 

Multiplying eg. (A. 11) by density and Favre-averaging one has: 


pH 


■P + pQ + p^- 


(A. 12) 


<P'P> and input of turbulent kinetic energy are 


v(r - iy 

where correlation 
neglected. 

The values of r/(r-l) and Q depend on reactive species mass fractions 
and temperature only. So, they can be calculated using flamelet solutions 
(A.7). To obtain their mean values, the averaging procedure of item A.4 was 


used. The value of H was obtained from the Favre-averaged energy 
conservation equation in its conventional form. As the result, eq.(A. 12) gave 
the relation between mean values of density p, pressure P and mean flow 

velocity U . Here, only two parameters (r and Q) were needed from the 
flamelet library to close hydrodynamics equations. The particular solutions 
for Q and Y at required valued of mixture fraction, scalar dissipation and 
pressure were obtained from the flamelet library by the linear interpolation 
between neighboring solutions in the same manner as it was done in Bezgin 
et.al. (1996). 
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A ppendix B Inflow Boundary Layers Simulation 


The inflows, in the considered test cases, had boundary layers. So the 
inflow distributions for CFD were modeled taking into account this feature. 
The longitudinal velocity and eddy viscosity distributions in the incoming 
flows were derived from the experimental data on boundary layer thickness 
5. The eddy viscosity distribution inside the boundary layer was 
approximated using modified Van-Driest relation: 

yU, y 

v, =0.41yU t (l-e 26v )(1 - ^)- e' s 

5 


where dynamics velocity U T was calculated using correlation for local skin 
friction factor c f for compressible turbulent boundary layer at flat plate: 


2U: 


U: 


= c, = 0.023 


U§ 


02 / 


(k-0 

+ 0.7^ -M' 


- 0.5 


— \ 0.5 


( l+T .) 


Here subscript "e" denotes parameters in core flow, T w is the temperature 
factor, k is heat capacities ratio. The longitudinal velocity distribution in the 
boundary layer was obtained by integrating of the equation for the shear 
stress: 


_ da 

=P(v t + v)— ; 

dy 

where the distribution for the shear stress inside the boundary layer was 
approximated as: 


( 


= pU; 


1-3 


+2i y 


-A 


The turbulence kinetic energy distribution in the boundary layer was 
ed from the approximate "equilibrium turbulence” relation: 

0.3; 


V, 

0U 

K 

dy 


The same procedure was used for u , v t and K distributions 
calculations (if required) in the exit plane of the hydrogen injectors. 
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A ppendix C Thermodynamics and Kinetics Reference Data 


Table C.l. 


Reaction rates constants for hydrogen oxidation chemistry based on data by 

Miller and Bowman (1989). 


n C 

Reaction rate constant is presented in the form K = AT exp(- . 

K 1 

Units used are Kelvins, seconds, mols and centimeters. 



Reaction 


Forward 

rate 

Ig A 

B 

-E/R 

1 

H2 + 02 = OH + OH 

13.23 

0.00 

-24044.0 

2 

OH + H2 = H20 + H 

9.07 

1.30 

-1824.7 

3 

O + OH = 02 + H 

14.60 

-0.50 

0.0 

4 

0 + H2 = OH + H 

4.70 

2.67 

-3165.3 

5 

H + 02 + M = H02 + M 

17.56 

-0.72 

0.0 

6 

OH + H02 = H20 + OH 

12.88 

0.00 

0.0 

7 

H + H02 = OH + OH 

14.15 

0.00 

-540.0 

8 

0 + H02 = 02 + OH 

13.15 

0.00 

-540.0 

9 

OH + OH = 0 + H20 

8.78 

1.30 

0.0 

10 

H + H + M = H2 + M 

18.00 

-1.00 

0.0 

11 

H + H + H2 = H2 + H2 

16.96 

-0.60 

0.0 

12 

H + H + H20 = H2 + H20 

19.78 

-1.25 

0.0 

13 

H + OH + M = H20 + M 

22.20 

-2.00 

0.0 

14 

H + O + M = OH + M 

16.79 

-0.60 

0.0 

15 

0 + 0 + M = 02 + M 

13.28 

0.00 

899.8 

16 

H + H02 = H2 + 02 

13.10 

0.00 

0.0 

17 

H02 + H02 = H202 + 02 

12.30 

0.00 

0.0 

18 

H202 + M = OH + OH + M 

17.11 

0.00 

-22896.7 

19 

H202 + H = H02 + H2 

12.20 

0.00 

-1912.2 

20 

H202 + OH = H20 + H02 

13.00 

0.00 

-905.8 

21 

H02 + NO = N02 + OH 

12.32 

0.00 

241.0 

22 

N02 + H = NO + OH 

14.54 

0.00 

-754.8 

23 

N02 + 0 = NO + 02 

13.00 

0.00 

-301.9 

24 

N02 + M = NO + O + M 

16.04 

0.00 

-33212.7 

25 

HNO + M = H + NO + M 

16.18 

0.00 

-24496.9 

26 

HNO + OH = NO + H20 

13.56 

0.00 

0.0 

27 

HNO + H = H2 + NO 

12.70 

0.00 

0.0 

28 

N + NO = N2 + 0 

12.51 

0.30 

0.0 

29 

N + 02 = NO + 0 

9.81 

1.00 

-3160.2 

30 

N + OH = NO + H 

13.58 

0.00 

0.0 


The third-body efficiencies are as follows: for reaction (5), H20 = 18.6, H2 =2.6 and N2 = 1.3; for 
reaction (10), H2 = 0.0 and H20 = 0.0; for reactions (13) and (14), H20 = 5.0; and for reaction (25), H20 
= 10.0, 02 = 2.0, N2 = 2.0 and H2 = 2 . 0 . 
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Table C.2 


Reaction rates constants for hydrogen oxidation chemistry based on data by 

Jachimowski (1988) 


n n 

Reaction rate constant is presented in the form K = AT exp(- — — ). 

R 1 


Units used are Kelvins, seconds, moles and centimeters. 



Reaction 


Forward 

rate 

IgA 

B 

-E/R 

1 

H2 + 02 = OH + OH 

13.23 

0.00 

-24154.7 

2 

H + 02 = OH + O 

14.41 

0.00 

-8454.1 

3 

0 + H2 = OH + H 

10.26 

1.00 

-4478.7 

4 

OH + H2 = H20 + H 

13.34 

0.00 

-2591.6 

5 

OH + OH = H20 + 0 

12.80 

0.00 

-548.5 

6 

H + OH + M = H20 + M 

22.34 

-2.00 

0.0 

7 

H + H + M = H2 + M 

17.81 

-1.00 

0.0 

8 

H + 0 + M = OH + M 

16.78 

-0.60 

0.0 

9 

H + 02 + M = H02 +M 

15.32 

0.00 

503.2 

10 

H02 + H = H2 + 02 

13.11 

0.00 

0.0 

11 

H02 + H = OH + OH 

14.15 

0.00 

-543.5 

12 

H02 + H = H20 + O 

13.00 

0.00 

-543.5 

13 

H02 + 0 = 02 + OH 

13.18 

0.00 

-478.1 

14 

H02 + OH = H20 + 02 

12.90 

0.00 

0.0 

15 

H02 + H02 = H202 + 02 

12.30 

0.00 

0.0 

16 

H + H202 = H2 + H02 

12.15 

0.00 

-1811.6 

17 

0 + H202 = OH + H02 

13.15 

0.00 

-4227.1 

18 

OH + H202 = H20 + H02 

12.79 

0.00 

-719.6 

19 

H202 + M = OH + OH +M 

17.08 

0.00 

-22896.7 

20 

0 + 0 + M = 02 + M 

17.78 

0.00 

905.8 

21 

N + N + M = N2 + M 

17.45 

-0.75 

0.0 

22 

N + 02 = NO + O 

9.81 

1.00 

-4176.8 

23 

N + NO = N2 + O 

13.20 

0.00 

0.0 

24 

N + OH = NO + H 

11.80 

0.50 

0.0 

25 

H + NO + M = HNO + M 

15.73 

0.00 

301.9 

26 

H + HNO = NO + H2 

12.68 

0.00 

0.0 

27 

0 + HNO = NO + OH 

11.70 

0.50 

0.0 

28 

OH + HNO = NO + H20 

13.56 

0.00 

0.0 

29 

H02 + HNO = NO + H202 

12.30 

0.00 

0.0 

30 

H02 + NO = N02 + OH 

12.53 

0.00 

130.8 

31 

H + N02 = NO + OH 

14.54 

0.00 

-754.8 

32 

O + N02 = NO + 02 

13.00 

0.00 

-301.9 

33 

N02 + M = NO + 0 + M 

16.06 

0.00 

-33212 7 


The third-body efficiencies relative to N2 = 1.0 are as follows: for reaction (6), H20 = 6.0; for reaction 
(7), H2 = 2.0 and H20 - 6.0; for reaction (8), H20 = 5.0; for reaction (9). H2 = 2.0 and H20 = 16.0; and 
reaction (19), H20= 15.0. 
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Polynomial coefficients for thermodynamical properties approximation from Alemasov, et.al. (1971). 
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/ FUEL LEAN SIDE (Z<Z g ) 



Fig. 1 . Schematic of the flamelet approach. 








FOUR PARAMETRIC DEPENDENCE Cot=(z,N,P,U 2 /2, boundary conditions) 






(H2/air, P = 0.05MPa, Tf = 900K, Ta = 900K) 



Fig. 3. Ignition/extinction behaviour of flamelet model. 



a) Beach experiment b) Burrows-Kurkov experiment 
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Fig. 4. Accuracy of previous FL predictions. 
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Fig. 5. Relative role of the chemical noncquilibriumness/mixture fraction fluctuations 
in prev ious FL approach predictions budget. 

Flamelel averaged solution. 

Flamelel instantaneous solution. 

Equilibrium chemistry limit 
^ Experimental data 

Burrows-Kurkov test case. 
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Fig. 6 Tests of detailed chemistry model. 

Miller & Bowman scheme (ignoring third body efficiences) 
Miller & Bowman scheme (selective third body efficiences) 
Jachimowski scheme 

Conditions of Burrows-Kurkov test case. 
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Fig. 7. Scheme and conditions of Burrows and Kurkov experiment 

*) Expected values are italic. They were obtained from equilibrium calculations. 
Species with mass fraction less than 10' 5 are not shown. 
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Fig. 9. Inflow parameters 
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Burrous-Kurkov test case 
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Fig. 1 1. Comparison of FL predictions with experimental data. 

x=0.356 m. 

current predictions 

— — - previous results 

+ experimental data 
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Burrows-Kurkov test case. 
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Fig. 14. Mixture fraction fluctuations intensity v.s. mean mixture fraction z 

x=0.356 m cross section. 

current modeling 

— — - previous results 


Burrows-Kurkov test case. 



Cut from computational region 


a) Mean mixture fraction: 



old calculations 
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b) Mixture fraction variance: 
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c) Turbulent kinetic energy: 
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Fig. 15. Turbulence characteristics 


Burrous-Kurkov test case 
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Fig. 16. Wall pressure distributions 

^ experimental data 

FL approach with small smoothing 
FL approach with large smoothing 
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Burrow s-Kurkov lest case 




Fig. 17 . Scheme and conditions of Beach et al. coaxial experiment. 

*) Expected values are italic . They were obtained from equilibrium calculations. 
Species with mass fraction less than 1(K 5 are not shown. 



Fig. 18. Schematic of computational domain 


Beach test case 
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Fig. 19. L2 norm behavior. 
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Beach test case 



Beach test case 
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Mach number contours 



Pressure contours 


Fig.22. Flow fields near injector 
(Domain I in Fig. 1 8) 


Beach test case 
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Fig. 25c. Comparison of FL predictions with experimental data. 

current predictions 

— — ■ previous results 

experimental data of Beach et al. 
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D = 0.01778 m 

Injector internal diameter d = 0.00236 m 
Injector lip thickness = 0.000725 m 
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Fig. 26. Scheme and conditions of Cheng et al. coaxial experiment. 
*) Expected values are italic. 



Fig. 27. Schematic of computational domain 


Cheng et.al. test case 
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Fig. 28. Inflow parameters distributions 


Cheng et.al. test case 
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Fig. 45. Profiles of turbulent kinetic energy. Non-reacting case. 
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Fig. 50. Profiles of turbulent kinetic energy. Reacting case. 
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